{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 18 - 异质干预效应与个性化\n",
    "\n",
    "## 从预测到因果推理\n",
    " \n",
    "在上一章中，我们简要介绍了机器学习模型。 ML 模型是用于我所说的预测的工具，或者更专业地说，是估计条件期望函数 \\\\(E[Y|X]\\\\)。换句话说，当您想从已知输入 \\\\(X\\\\) （如英语句子、本月销售、大脑扫描图像）映射到最初未知但定义明确的输出 \\\\(Y \\\\)（如日语句子、下个月的销售额或癌症诊断）。因此，如果 ML 处理预测或估计 \\\\(E[Y|X]\\\\)，为了使其有用，您必须将您想用 ML 解决的任何问题框定为预测问题，即估计 \\\\(E[Y|X]\\\\) 是关键。我们在上一章中介绍了这样一个例子。在那里，我们必须根据客户的特定特征来预测客户的盈利能力：\\\\(E[NetValue|Age, Income, Region]\\\\)。这些信息非常有用，因为它使我们能够集中精力与盈利客户打交道，而不是与非盈利客户开展业务。在这里，很好地预测盈利能力是关键。\n",
    " \n",
    "请注意，在您将自己从数据生成过程中移除的意义上，这是一种被动的估计方法。在我们的示例中，我们假设给出了客户盈利能力“净值”。我们所要做的就是估计它。换句话说，我们假设除了预测客户的盈利能力外，我们无能为力。我们不能增加它，也不能减少它。但这并不总是正确的。事实上，很多时候，公司都有可以用来提高客户盈利能力的杠杆。这些杠杆的范围可以从优质或更便宜的客户服务到折扣、价格或营销。在行业中，经常会出现我们被插入到数据生成过程中的情况。我们可以影响它。因此，作为在该行业工作的数据科学家，我们经常必须回答最佳行动方案或干预措施，以优化某些业务指标，通常是盈利能力或其他一些中间指标，如转换、成本或销售额。\n",
    " \n",
    "在这个我们不是被动观察的世界中，估计 \\\\(E[Y|X]\\\\) 并不是全部。这是我们进入因果推理的地方。我们需要在条件期望函数中添加另一部分，也就是对我们参与数据生成过程进行建模的部分，即干预本身：\n",
    " \n",
    "$$\n",
    "E[Y|X, T]\n",
    "$$\n",
    " \n",
    "我们现在必须区分上下文或外生特征\\\\(X\\\\) 和处理\\\\(T\\\\)。两者都会影响结果 \\\\(Y\\\\)，但是虽然我们无法控制 \\\\(X\\\\)，但我们可以决定 \\\\(T\\\\) 将取什么值，或者至少对其进行干预。举个具体的例子，\\\\(Y\\\\) 可能是一天的销售额，\\\\(X\\\\) 可能是您无法控制的上下文特征，但它会为您提供有关销售的信息，例如前几天，\\\\(T\\\\) 是您可以干预以增加销售额的处理变量，例如价格、商品库存水平或营销。因果推理是在上下文 \\\\(X\\\\) 下估计 \\\\(T\\\\) 和 \\\\(Y\\\\) 之间因果关系的过程。一旦我们这样做了，优化 \\\\(Y\\\\) 只是以最佳方式设置处理 \\\\(T\\\\) 的问题\n",
    " \n",
    "$$\n",
    "\\underset{T}{argmax} \\ E[Y|X, T]\n",
    "$$\n",
    " \n",
    " \n",
    "从这个意义上说，除了因果推理的积极方面，我们还有一个规范的动机。\n",
    " \n",
    "在第一部分，我们试图回答诸如学校教育的价值是什么？法律变化可以降低吸烟水平吗？我们可以通过积极的心态来提高学业成绩吗？酒精对死亡率的影响是什么？从理解世界如何运作的纯科学观点来看，所有这些问题都很有趣。但它们背后也有实际的动机。如果我们知道学校教育对收入的影响，我们就可以理解为此付出的合理价格。用数学术语来说，我们所做的是估计学校教育的因果推理并对其进行优化：\\\\(\\underset{Educ}{argmax} \\ E[Income|X, Educ]\\\\)。\n",
    " \n",
    "第一部分的重点是回答干预总体上是积极的、强的还是零。例如，我们想知道一般而言，投资于教育是否是一个好主意。同样在第一部分中，\\\\(X\\\\) 的作用是双重的。首先，\\\\(X\\\\) 可能包含混杂因素，在这种情况下，因果效应只有在我们考虑或调整 \\\\(X\\\\) 时才能识别。或者，\\\\(X\\\\) 可以减少因果估计的方差。如果 \\\\(X\\\\) 是 \\\\(Y\\\\) 的良好预测器，我们可以用它来解释 \\\\(Y\\\\) 的方差，从而使因果效应更加明显。\n",
    " \n",
    "现在，事情将变得不那么黑白分明了。我们想要的不仅仅是平均干预效果。我们将允许干预对某些人产生积极影响，但对其他人则不然。上下文特征 \\\\(X\\\\) 将在定义不同的单位配置文件中发挥作用，每个配置文件可能对处理的反应不同。我们现在想要"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 从 ATE 到 CATE\n",
    " \n",
    " \n",
    "到目前为止，我们每次估计干预的因果影响时，都是平均干预效果（或者有时是局部平均干预效果）：\n",
    " \n",
    "$$\n",
    "E[Y_1−Y_0]\n",
    "$$\n",
    " \n",
    "或等效的连续干预\n",
    " \n",
    "$$\n",
    "E[y'(t)]\n",
    "$$\n",
    " \n",
    "其中 \\\\(y'(t)\\\\) 是响应函数或结果的处理导数。我们已经学会了揭示干预的一般有效性的技术。 ATE 估计是因果推理的基础。对于我们称为程序评估的决策问题，它是一个超级有用的工具。我们想知道我们是否应该向整个人群推出一种干预方法。不要被公共政策条款所迷惑。评估国家教育或健康计划有效性的相同技术也可用于了解推出新产品对公司底线的影响。这里要注意的关键是我们想要告知的决定是我们是否应该干预。\n",
    " \n",
    "现在，我们将尝试告知另一种类型的决定：**我们干预谁**？现在，我们允许决定从一个单元更改为另一个单元。干预一个单位而不是另一个单位可能是有益的。我们希望个性化干预。用更专业的术语来说，我们想要估计条件平均干预效果 (CATE)\n",
    " \n",
    "$$\n",
    "E[Y_1−Y_0 | X] \\ \\text{或} \\ E[y'(t)|X]\n",
    "$$\n",
    " \n",
    "对 \\\\(X\\\\) 的条件化意味着我们现在允许处理效果根据每个单元的特性而有所不同。同样，在这里，我们认为并非所有实体对干预的反应都一样好。我们希望利用这种异质性。我们只想处理正确的单位（在二元情况下）或弄清楚每个单位的最佳干预程度是多少（在连续情况下）。\n",
    " \n",
    "例如，如果您是一家必须决定每个客户有资格获得贷款的银行，那么您可以确定向所有人提供大量资金并不是一个好主意——尽管这对某些人来说可能是合理的。您必须明智地对待您的干预手段（贷款金额）。也许，根据客户的信用评分 (\\\\(X\\\\))，您可以找出合适的贷款剂量。当然，您无需成为大型机构即可利用个性化。不乏适用的例子。您应该在一年中的哪几天进行销售？你应该为任何产品收取多少费用？对每个人来说，多少运动量才算过量运动量？\n",
    " \n",
    "这样想吧。您有一群客户和一种待遇（价格、折扣、贷款……）。您想要个性化的待遇，例如，给不同的客户不同的折扣。\n",
    " \n",
    "![img](./data/img/causal-model/customers.png)\n",
    " \n",
    "为此，您必须对客户进行细分。您创建了对您的干预有不同反应的小组。例如，您希望找到对折扣反应良好的客户和对折扣反应不佳的客户。好吧，客户对处理的反应由条件处理效果 \\\\(\\frac{\\delta Y}{ \\delta T}\\\\) 给出。因此，我们可以以某种方式估计，对于每个客户，我们可以将那些对干预反应很好（高干预效果）和那些对干预反应不佳的人分组在一起。如果我们这样做，我们会像下图那样拆分客户空间。\n",
    " \n",
    "![img](./data/img/causal-model/elast-partition.png)\n",
    " \n",
    "这将是美妙的，因为现在我们将能够估计每个分区上的不同处理效果或弹性。请注意，弹性只是从 \\\\(T\\\\) 到 \\\\(Y\\\\) 的直线或函数的斜率。因此，如果我们可以生成斜率或弹性不同的分区，这意味着这些分区上的实体对处理具有不同的响应性。\n",
    " \n",
    "![img](./data/img/causal-model/elast-split.png)\n",
    " \n",
    "换句话说，您想要摆脱以原始形式预测 \\\\(Y\\\\) 并开始预测 \\\\(Y\\\\) 在 \\\\(T\\\\) 上的导数，\\\\( \\frac{\\delta Y}{ \\delta T}\\\\) 为每个单位。例如，假设 \\\\(Y\\\\) 是冰淇淋销量，\\\\(T\\\\) 是冰淇淋价格，每个单位 \\\\(i\\\\) 是一天。让我们把道德问题放在一边，为了争论，假装你每天都能改变冰淇淋的价格。如果您能以某种方式找到 \\\\(\\frac{\\delta Sales}{ \\delta Price}\\\\) **低的日子**，那么您可以**提高价格** 而不会在那些日子损失太多销售额。也许你已经这样做了，比如说，当你在假期增加它们时。关键是，根据价格弹性来区分日期是很有用的，因为它为您提供了如何以最佳方式设定价格的基础。\n",
    " \n",
    "好吧，你可能会说，但这有点棘手。如果我看不到弹性 \\\\(\\frac{\\delta Sales}{ \\delta Price}\\\\)，如何预测它？这是一个很好的观点。弹性在单位级别上基本上是不可观察的。不仅如此，这是一个奇怪的概念。我们更习惯于从原始数量的角度来思考，而不是根据这些相同数量的变化率来思考。因此，为了更好地概念化弹性，这里有一个小技巧。您可以将每个实体视为具有 \\\\(Y_i\\\\) 值，在我们的示例中为销售额，但也具有个体弹性 \\\\(\\frac{delta Y_i}{\\delta T_i}\\\\)。弹性是 \\\\(Y\\\\) 随 \\\\(T\\\\) 变化的程度，因此您可以考虑每个实体也具有与其相关联的斜率系数 \\\\(\\frac{\\delta Y}{ \\delta T}_i\\\\)。在我们的示例中，我们会说每一天的销售价格都有一个斜率系数。\n",
    " \n",
    "![img](./data/img/causal-model/elasticity.png)\n",
    " \n",
    "当然，我们看不到那些单独的斜率系数。为了让我们看到各个斜率，我们必须每天在两个不同的价格下观察，并计算每个价格的销售额如何变化。\n",
    " \n",
    "$$\n",
    "\\frac{\\delta Y_i}{ \\delta T_i} \\近似 \\frac{Y(T_i) - Y(T_i + \\epsilon)}{T_i - (T_i + \\epsilon)}\n",
    "$$\n",
    " \n",
    "这又是因果推理的根本问题。我们永远无法在不同的处理条件下看到相同的单元。所以，我们能做些什么？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 预测弹性\n",
    " \n",
    "我们在这里陷入了复杂的境地。我们已经同意我们需要预测 \\\\(\\frac{\\delta Y_i}{ \\delta T_i}\\\\)，遗憾的是这是不可观察的。因此，我们不能使用 ML 算法并将其作为目标插入。但也许我们不需要观察 \\\\(\\frac{\\delta Y_i}{ \\delta T_i}\\\\) 来预测它\n",
    " \n",
    "这是一个想法。如果我们使用线性回归呢？\n",
    " \n",
    "![img](./data/img/causal-model/linear-fix.png)\n",
    " \n",
    "假设您将以下线性模型拟合到您的数据中。\n",
    " \n",
    "$$\n",
    "y_i = \\beta_0 + \\beta_1 t_i + \\beta_2 X_i + e_i\n",
    "$$\n",
    " \n",
    "如果你在干预变量上区分它，你最终会得到\n",
    " \n",
    "$$\n",
    "\\frac{\\delta y_i}{\\delta t_i} = \\beta_1\n",
    "$$\n",
    " \n",
    "既然你可以估计上面的模型得到 \\\\(\\hat{\\beta_1}\\\\)，我们甚至可以大胆地说**即使你无法观察到弹性，你也可以预测它**。在上面的例子中，这是一个相当简单的预测，即我们为每个人预测常数值 \\\\(\\hat{\\beta_1}\\\\)。那是一些东西，但还不是我们想要的。那是ATE，不是CATE。这无助于我们根据实体对干预的反应程度对实体进行分组的任务，因为每个人都得到相同的弹性预测。但是，我们可以做以下简单的改变\n",
    " \n",
    "$$\n",
    "y_i = \\beta_0 + \\beta_1 t_i + \\beta_2 X_i + \\beta_3 t_i X_i + e_i\n",
    "$$\n",
    " \n",
    "这反过来会给我们以下弹性预测\n",
    " \n",
    "$$\n",
    "\\widehat{\\frac{\\delta y_i}{\\delta t_i}} = \\hat{\\beta_1} + \\hat{\\beta_3}X_i\n",
    "$$\n",
    " \n",
    "其中 \\\\(\\beta_3\\\\) 是 \\\\(X\\\\) 中特征的向量系数。\n",
    " \n",
    "现在由不同的 \\\\(X_i\\\\) 定义的每个实体将有不同的弹性预测。换句话说，弹性预测会随着 \\\\(X\\\\) 的变化而变化。唉，回归可以为我们提供一种估计 CATE \\\\(E[y'(t)|X]\\\\) 的方法。\n",
    " \n",
    "我们终于到了某个地方。上面的模型允许我们对每个实体进行弹性预测。通过这些预测，我们可以创建更多有用的组。我们可以将具有高预测弹性的单元组合在一起。我们可以对预测弹性低的那些做同样的事情。最后，通过我们的弹性预测，我们可以根据我们认为实体对干预的反应程度对它们进行分组。\n",
    " \n",
    "理论到此为止。是时候通过一个示例来说明如何制作这种弹性模型了。让我们考虑一下我们的冰淇淋示例。每个单位 \\\\(i\\\\) 是一天。对于每一天，我们都知道是否是工作日，制作冰淇淋的成本是多少（您可以将成本视为质量的代表）以及当天的平均温度。这些将是我们的特征空间 \\\\(X\\\\)。然后，我们有我们的干预手段、价格和我们的结果，即售出的冰淇淋数量。对于这个例子，我们将考虑随机处理，这样我们现在不必担心偏差。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import seaborn as sns\n",
    "import statsmodels.formula.api as smf\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5000, 5)\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>temp</th>\n",
       "      <th>weekday</th>\n",
       "      <th>cost</th>\n",
       "      <th>price</th>\n",
       "      <th>sales</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>25.8</td>\n",
       "      <td>1</td>\n",
       "      <td>0.3</td>\n",
       "      <td>7</td>\n",
       "      <td>230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>22.7</td>\n",
       "      <td>3</td>\n",
       "      <td>0.5</td>\n",
       "      <td>4</td>\n",
       "      <td>190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>33.7</td>\n",
       "      <td>7</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5</td>\n",
       "      <td>237</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>23.0</td>\n",
       "      <td>4</td>\n",
       "      <td>0.5</td>\n",
       "      <td>5</td>\n",
       "      <td>193</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>24.4</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>252</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   temp  weekday  cost  price  sales\n",
       "0  25.8        1   0.3      7    230\n",
       "1  22.7        3   0.5      4    190\n",
       "2  33.7        7   1.0      5    237\n",
       "3  23.0        4   0.5      5    193\n",
       "4  24.4        1   1.0      3    252"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prices_rnd = pd.read_csv(\"./data/ice_cream_sales_rnd.csv\")\n",
    "print(prices_rnd.shape)\n",
    "prices_rnd.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "请记住我们的目标：我们需要根据具体的“功能”、“临时”、“工作日”和“成本”来决定何时收取更多费用以及何时收取更少费用。 如果这是目标，则需要评估干预效果异质性模型在实现该目标方面的有用性。 我们稍后会谈到这一点（下一章会详细介绍）。 现在，让我们将数据集拆分为训练和测试集。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "train, test = train_test_split(prices_rnd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们有了训练数据，我们需要建立一个模型，让我们能够区分价格弹性高的日子和价格弹性低的日子。 我们这样做的方法是简单地预测价格弹性。 具体如何？ 首先，让我们考虑使用以下线性模型\n",
    " \n",
    "$$\n",
    "sales_i = \\beta_0 + \\beta_1 price_i + \\pmb{\\beta_2}X_i + e_i\n",
    "$$\n",
    " \n",
    "如果我们检查这个模型的参数，我们可以看到我们预测的弹性会是什么样子。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>  186.7113</td> <td>    1.770</td> <td>  105.499</td> <td> 0.000</td> <td>  183.241</td> <td>  190.181</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.2]</th> <td>  -25.0512</td> <td>    0.924</td> <td>  -27.114</td> <td> 0.000</td> <td>  -26.863</td> <td>  -23.240</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.3]</th> <td>  -24.5834</td> <td>    0.901</td> <td>  -27.282</td> <td> 0.000</td> <td>  -26.350</td> <td>  -22.817</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.4]</th> <td>  -24.3807</td> <td>    0.897</td> <td>  -27.195</td> <td> 0.000</td> <td>  -26.138</td> <td>  -22.623</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.5]</th> <td>  -24.9036</td> <td>    0.894</td> <td>  -27.850</td> <td> 0.000</td> <td>  -26.657</td> <td>  -23.150</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.6]</th> <td>  -24.0921</td> <td>    0.903</td> <td>  -26.693</td> <td> 0.000</td> <td>  -25.862</td> <td>  -22.323</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.7]</th> <td>   -0.8635</td> <td>    0.888</td> <td>   -0.972</td> <td> 0.331</td> <td>   -2.605</td> <td>    0.878</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>price</th>           <td>   -2.7515</td> <td>    0.106</td> <td>  -25.970</td> <td> 0.000</td> <td>   -2.959</td> <td>   -2.544</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>temp</th>            <td>    1.9848</td> <td>    0.060</td> <td>   33.117</td> <td> 0.000</td> <td>    1.867</td> <td>    2.102</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cost</th>            <td>    4.4718</td> <td>    0.528</td> <td>    8.462</td> <td> 0.000</td> <td>    3.436</td> <td>    5.508</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m1 = smf.ols(\"sales ~ price + temp+C(weekday)+cost\", data=train).fit()\n",
    "m1.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于 \\\\(m1\\\\)，预测的价格弹性 \\\\(\\widehat{\\dfrac{\\delta y_i}{\\delta t_i}}\\\\) 将由 \\\\(\\hat{\\beta_1}\\\\ )，在我们的例子中是 -2.75。这意味着，我们为冰淇淋收取的每额外巴西雷亚尔，我们预计销售额将下降约 3 个单位。\n",
    " \n",
    "注意这个 \\\\(m1\\\\) 如何为每个人预测完全相同的弹性。因此，如果我们想知道人们在哪些日子对冰淇淋价格不那么敏感，这不是一个很好的模型。当我们在这里需要的是 CATE 时，它会估计 ATE。请记住，我们的目标是以这样一种方式对实体进行分区，以便我们可以针对每个单独的分区个性化和优化我们的处理（价格）。如果每个预测都相同，我们就无法进行分区。我们没有区分敏感单位和非敏感单位。为了纠正这一点，考虑我们的第二个模型：\n",
    " \n",
    "$$\n",
    "sales_i = \\beta_0 + \\beta_1 price_i + \\beta_2 price_i * temp_i * + \\pmb{\\beta_3}X_i + e_i\n",
    "$$\n",
    " \n",
    "第二个模型包括价格和温度之间的**交互项**。这意味着它允许弹性在不同温度下有所不同。我们在这里实际上要说的是，人们对价格上涨或多或少敏感，具体取决于温度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "         <td></td>            <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>       <td>  192.4767</td> <td>    4.371</td> <td>   44.037</td> <td> 0.000</td> <td>  183.907</td> <td>  201.046</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.2]</th> <td>  -25.0805</td> <td>    0.924</td> <td>  -27.143</td> <td> 0.000</td> <td>  -26.892</td> <td>  -23.269</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.3]</th> <td>  -24.5871</td> <td>    0.901</td> <td>  -27.290</td> <td> 0.000</td> <td>  -26.354</td> <td>  -22.821</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.4]</th> <td>  -24.4225</td> <td>    0.897</td> <td>  -27.231</td> <td> 0.000</td> <td>  -26.181</td> <td>  -22.664</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.5]</th> <td>  -24.8953</td> <td>    0.894</td> <td>  -27.844</td> <td> 0.000</td> <td>  -26.648</td> <td>  -23.142</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.6]</th> <td>  -24.1269</td> <td>    0.903</td> <td>  -26.726</td> <td> 0.000</td> <td>  -25.897</td> <td>  -22.357</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>C(weekday)[T.7]</th> <td>   -0.8581</td> <td>    0.888</td> <td>   -0.966</td> <td> 0.334</td> <td>   -2.599</td> <td>    0.883</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>price</th>           <td>   -3.6299</td> <td>    0.618</td> <td>   -5.873</td> <td> 0.000</td> <td>   -4.842</td> <td>   -2.418</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>temp</th>            <td>    1.7459</td> <td>    0.176</td> <td>    9.912</td> <td> 0.000</td> <td>    1.401</td> <td>    2.091</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>price:temp</th>      <td>    0.0366</td> <td>    0.025</td> <td>    1.443</td> <td> 0.149</td> <td>   -0.013</td> <td>    0.086</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>cost</th>            <td>    4.4558</td> <td>    0.529</td> <td>    8.431</td> <td> 0.000</td> <td>    3.420</td> <td>    5.492</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.table.SimpleTable'>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m2 = smf.ols(\"sales ~ price*temp + C(weekday) + cost\", data=train).fit()\n",
    "m2.summary().tables[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一旦我们估计模型，预测的弹性由下式给出\n",
    " \n",
    "$$\n",
    "\\widehat{\\frac{\\delta sales_i}{\\delta price_i}} = \\hat{\\beta_1} + \\hat{\\beta_3}temp_i\n",
    "$$\n",
    " \n",
    "请注意，\\\\(\\hat{\\beta_3}\\\\) 为正 0,03，基线弹性 \\\\(\\beta_1\\\\)（在 \\\\(0C^o\\\\) 处的弹性）为 -3.6。这意味着，平均而言，随着我们提高价格，销售额下降，这是有道理的。这也意味着温度每升高一度，人们对冰淇淋价格上涨的敏感度就会降低（尽管幅度不大）。例如，在 \\\\(25C^o\\\\)，我们每多收取一次巴西雷亚尔，我们的销售额就会下降 2.8 个单位 \\\\((-3.6 + (0.03 * 25))\\\\)。但是在 \\\\(35C^o\\\\)，我们每多收取一次巴西雷亚尔，它们只会下降 2.5 个单位 \\\\((-3.6 + (0.03 * 35))\\\\)。这也是一种直觉。随着天气越来越热，人们愿意为冰淇淋支付更多费用。\n",
    " \n",
    "我们可以走得更远。下一个模型包括所有特征空间的交互项。这意味着弹性会随着温度、星期几和成本而变化。\n",
    " \n",
    "$$\n",
    "sales_i = \\beta_0 + \\beta_1 price_i + \\pmb{\\beta_2 X_i}*price_i + \\pmb{\\beta_3}X_i + e_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "m3 = smf.ols(\"sales ~ price*cost + price*C(weekday) + price*temp\", data=train).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据上述模型，单位水平弹性或 CATE 将由下式给出\n",
    " \n",
    "$$\n",
    "\\frac{\\delta 销售额}{\\delta 价格} = \\beta_1 + \\pmb{\\beta_2 X_i}\n",
    "$$\n",
    " \n",
    "其中 \\\\(\\beta_1\\\\) 是价格系数，\\\\(\\pmb{\\beta_2}\\\\) 是交互系数的向量。\n",
    " \n",
    "最后，让我们看看如何实际做出这些弹性预测。一种方法是从模型中提取弹性参数并使用上面的公式。但是，我们将采用更一般的近似值。由于弹性只不过是干预结果的导数，我们可以使用导数的定义。\n",
    " \n",
    "$$\n",
    "\\frac{\\delta y}{\\delta t} = \\dfrac{y(t+\\epsilon) - y(t)}{ (t + \\epsilon) - t }\n",
    "$$\n",
    " \n",
    "随着 \\\\(\\epsilon\\\\) 变为零。我们可以通过将 \\\\(\\epsilon\\\\) 替换为 1 来近似这个定义。\n",
    " \n",
    "$$\n",
    "\\frac{\\delta y}{\\delta t} \\近似 \\hat{y}(t+1) - \\hat{y}(t)\n",
    "$$\n",
    " \n",
    "其中 \\\\(\\hat{y}\\\\) 由我们模型的预测给出。换句话说，我将用我的模型做出两个预测：一个传递原始数据，另一个传递原始数据，但处理增加一个单位。这些预测之间的差异是我的 CATE 预测。\n",
    " \n",
    "下面，您可以看到执行此操作的函数。由于我们使用了训练集来估计我们的模型，我们现在将对测试集进行预测。首先，让我们使用我们的第一个 ATE 模型 \\\\(m1\\\\)。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>temp</th>\n",
       "      <th>weekday</th>\n",
       "      <th>cost</th>\n",
       "      <th>price</th>\n",
       "      <th>sales</th>\n",
       "      <th>pred_elast</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2648</th>\n",
       "      <td>18.6</td>\n",
       "      <td>7</td>\n",
       "      <td>0.5</td>\n",
       "      <td>10</td>\n",
       "      <td>185</td>\n",
       "      <td>-2.751463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2456</th>\n",
       "      <td>26.0</td>\n",
       "      <td>3</td>\n",
       "      <td>0.5</td>\n",
       "      <td>10</td>\n",
       "      <td>200</td>\n",
       "      <td>-2.751463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4557</th>\n",
       "      <td>23.7</td>\n",
       "      <td>3</td>\n",
       "      <td>0.3</td>\n",
       "      <td>8</td>\n",
       "      <td>192</td>\n",
       "      <td>-2.751463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4884</th>\n",
       "      <td>28.9</td>\n",
       "      <td>4</td>\n",
       "      <td>1.5</td>\n",
       "      <td>6</td>\n",
       "      <td>213</td>\n",
       "      <td>-2.751463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>92</th>\n",
       "      <td>23.7</td>\n",
       "      <td>1</td>\n",
       "      <td>0.5</td>\n",
       "      <td>8</td>\n",
       "      <td>207</td>\n",
       "      <td>-2.751463</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      temp  weekday  cost  price  sales  pred_elast\n",
       "2648  18.6        7   0.5     10    185   -2.751463\n",
       "2456  26.0        3   0.5     10    200   -2.751463\n",
       "4557  23.7        3   0.3      8    192   -2.751463\n",
       "4884  28.9        4   1.5      6    213   -2.751463\n",
       "92    23.7        1   0.5      8    207   -2.751463"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pred_elasticity(m, df, t=\"price\"):\n",
    "    return df.assign(**{\n",
    "        \"pred_elast\": m.predict(df.assign(**{t:df[t]+1})) - m.predict(df)\n",
    "    })\n",
    "\n",
    "pred_elasticity(m1, test).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用 \\\\(m1\\\\) 进行弹性预测并不好玩。 我们可以看到它预测了所有日子的完全相同的值。 那是因为该模型上没有交互项。 但是，如果我们使用 \\\\(m3\\\\) 进行预测，它会为每天输出不同的弹性预测。 那是因为现在，弹性或干预效果取决于当天的具体特征。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>temp</th>\n",
       "      <th>weekday</th>\n",
       "      <th>cost</th>\n",
       "      <th>price</th>\n",
       "      <th>sales</th>\n",
       "      <th>pred_elast</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4764</th>\n",
       "      <td>31.1</td>\n",
       "      <td>6</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3</td>\n",
       "      <td>212</td>\n",
       "      <td>1.144309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4324</th>\n",
       "      <td>24.8</td>\n",
       "      <td>7</td>\n",
       "      <td>0.5</td>\n",
       "      <td>10</td>\n",
       "      <td>182</td>\n",
       "      <td>-9.994303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4536</th>\n",
       "      <td>25.0</td>\n",
       "      <td>2</td>\n",
       "      <td>1.5</td>\n",
       "      <td>6</td>\n",
       "      <td>205</td>\n",
       "      <td>0.279273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3466</th>\n",
       "      <td>26.0</td>\n",
       "      <td>3</td>\n",
       "      <td>1.5</td>\n",
       "      <td>3</td>\n",
       "      <td>205</td>\n",
       "      <td>0.308320</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>115</th>\n",
       "      <td>19.3</td>\n",
       "      <td>3</td>\n",
       "      <td>0.3</td>\n",
       "      <td>9</td>\n",
       "      <td>177</td>\n",
       "      <td>-0.349745</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      temp  weekday  cost  price  sales  pred_elast\n",
       "4764  31.1        6   1.0      3    212    1.144309\n",
       "4324  24.8        7   0.5     10    182   -9.994303\n",
       "4536  25.0        2   1.5      6    205    0.279273\n",
       "3466  26.0        3   1.5      3    205    0.308320\n",
       "115   19.3        3   0.3      9    177   -0.349745"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pred_elast3 = pred_elasticity(m3, test)\n",
    "\n",
    "np.random.seed(1)\n",
    "pred_elast3.sample(5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "请注意预测是如何从 -9 变为 1 的数字。这些不是销售列的预测，其数量级为数百。相反，**这是一个预测，如果我们将价格提高一个单位，销售额会发生多少变化**。在赌注之外，我们可以看到一些奇怪的数字。例如，看一下第 4764 天。它预测的是正弹性。换句话说，我们预测如果我们提高冰淇淋的价格，销售额将会增加。这不符合我们的经济意识。可能是模型对该预测进行了一些奇怪的外推。幸运的是，您不必为此担心太多。请记住，我们的最终目标是根据单位对干预的敏感程度来划分单位。 **不是**提出有史以来最准确的弹性预测。对于我们的主要目标，如果弹性预测根据单元的敏感程度对单元进行排序就足够了。换句话说，即使像 1.1 或 0.5 这样的正弹性预测没有多大意义，我们所需要的只是排序正确，也就是说，我们希望预测为 1.1 的单位比单位受价格上涨的影响更小预测为 0.5。\n",
    " \n",
    "好的，我们有弹性或 CATE 模型。但仍然存在一个潜在的问题：它们与 ML 预测模型相比如何？现在让我们试试吧。我们将使用机器学习算法，使用价格、温度、工作日和成本作为特征 \\\\(X\\\\) 并尝试预测冰淇淋销售。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9124088322890127"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = [\"temp\", \"weekday\", \"cost\", \"price\"]\n",
    "y = \"sales\"\n",
    "ml = GradientBoostingRegressor()\n",
    "ml.fit(train[X], train[y])\n",
    "\n",
    "# make sure the model is not overfiting.\n",
    "ml.score(test[X], test[y])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该模型可以预测我们每天的销售额。 但它适合我们真正想要的吗？ 换句话说，这个模型能区分人们对冰淇淋价格更敏感的日子吗？ 它可以帮助我们根据价格敏感度决定收取多少费用吗？\n",
    " \n",
    "要查看哪个模型更有用，让我们尝试使用它们来分割个体单元。 对于每个模型，我们将个体单元分成 2 组。 我们希望一个群体对价格上涨反应灵敏，而另一个群体则反应不大。 如果是这样的话，我们可以围绕这些群体组织我们的业务：对于属于高响应群体的日子，我们最好不要将价格定得太高。 对于低响应群体，我们可以提高价格而不会冒太大的销售风险。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>temp</th>\n",
       "      <th>weekday</th>\n",
       "      <th>cost</th>\n",
       "      <th>price</th>\n",
       "      <th>sales</th>\n",
       "      <th>pred_elast</th>\n",
       "      <th>elast_band</th>\n",
       "      <th>pred_sales</th>\n",
       "      <th>pred_band</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2648</th>\n",
       "      <td>18.6</td>\n",
       "      <td>7</td>\n",
       "      <td>0.5</td>\n",
       "      <td>10</td>\n",
       "      <td>185</td>\n",
       "      <td>-10.301045</td>\n",
       "      <td>(-10.597999999999999, -0.00555]</td>\n",
       "      <td>186.878081</td>\n",
       "      <td>(161.089, 198.735]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2456</th>\n",
       "      <td>26.0</td>\n",
       "      <td>3</td>\n",
       "      <td>0.5</td>\n",
       "      <td>10</td>\n",
       "      <td>200</td>\n",
       "      <td>0.036165</td>\n",
       "      <td>(-0.00555, 1.389]</td>\n",
       "      <td>203.188327</td>\n",
       "      <td>(198.735, 257.746]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4557</th>\n",
       "      <td>23.7</td>\n",
       "      <td>3</td>\n",
       "      <td>0.3</td>\n",
       "      <td>8</td>\n",
       "      <td>192</td>\n",
       "      <td>-0.132057</td>\n",
       "      <td>(-10.597999999999999, -0.00555]</td>\n",
       "      <td>188.800637</td>\n",
       "      <td>(161.089, 198.735]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4884</th>\n",
       "      <td>28.9</td>\n",
       "      <td>4</td>\n",
       "      <td>1.5</td>\n",
       "      <td>6</td>\n",
       "      <td>213</td>\n",
       "      <td>0.860663</td>\n",
       "      <td>(-0.00555, 1.389]</td>\n",
       "      <td>210.430813</td>\n",
       "      <td>(198.735, 257.746]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>92</th>\n",
       "      <td>23.7</td>\n",
       "      <td>1</td>\n",
       "      <td>0.5</td>\n",
       "      <td>8</td>\n",
       "      <td>207</td>\n",
       "      <td>-9.953698</td>\n",
       "      <td>(-10.597999999999999, -0.00555]</td>\n",
       "      <td>209.044522</td>\n",
       "      <td>(198.735, 257.746]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      temp  weekday  cost  price  sales  pred_elast  \\\n",
       "2648  18.6        7   0.5     10    185  -10.301045   \n",
       "2456  26.0        3   0.5     10    200    0.036165   \n",
       "4557  23.7        3   0.3      8    192   -0.132057   \n",
       "4884  28.9        4   1.5      6    213    0.860663   \n",
       "92    23.7        1   0.5      8    207   -9.953698   \n",
       "\n",
       "                           elast_band  pred_sales           pred_band  \n",
       "2648  (-10.597999999999999, -0.00555]  186.878081  (161.089, 198.735]  \n",
       "2456                (-0.00555, 1.389]  203.188327  (198.735, 257.746]  \n",
       "4557  (-10.597999999999999, -0.00555]  188.800637  (161.089, 198.735]  \n",
       "4884                (-0.00555, 1.389]  210.430813  (198.735, 257.746]  \n",
       "92    (-10.597999999999999, -0.00555]  209.044522  (198.735, 257.746]  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bands_df = pred_elast3.assign(\n",
    "    elast_band = pd.qcut(pred_elast3[\"pred_elast\"], 2), # create two groups based on elasticity predictions \n",
    "    pred_sales = ml.predict(pred_elast3[X]),\n",
    "    pred_band = pd.qcut(ml.predict(pred_elast3[X]), 2), # create two groups based on sales predictions\n",
    ")\n",
    "\n",
    "bands_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们需要比较这两个人群分层模型哪一个是最好的。 我现在可能有点超前了，因为我们只会在下一章中讨论 CATE 模型评估。 但我觉得我可以让你尝尝它的样子。 检查这些分层模式有多好的一种非常简单的方法 - 我的意思是好用 - 就是绘制每个分区的销售价格回归线。 我们可以通过 Seaborn 的 `regplot` 和 `FacetGrid` 轻松实现这一点。\n",
    " \n",
    "下面，我们可以看到使用弹性预测进行的分层。 请记住，所有这些都是在测试集中完成的。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAADQCAYAAABStPXYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABG60lEQVR4nO29eZxcVZnw/31q6707naWzdjZICIsBQkAQhIiouIE6jOL7jqOODujwCjqjo4wj4zCDLy6DA+M7YxjGbUZBxC36A9SgMUZZTEL2BBI6W2fr7qTTa3Vt9/n9cW51V3dXpbtS93ZXdZ/v51OfuvXce0+duz7nPOc5zyOqisVisVgsxUZgvCtgsVgsFks2rIKyWCwWS1FiFZTFYrFYihKroCwWi8VSlFgFZbFYLJaixCooi8VisRQlVkFZLBaLpSjxXEGJSEpEtmR8PuPK14nIyrMo7x0icsEotlslIh3uf24TkbUi0nA2x5Cl7A+IyNfOUL97cqy7T0QOi0j3EHmZiHxfRPaJyPMisjDH/utE5KWMc9ngyheIyDPuca4TkXmu/HVDzn2fiLzDXXe9iGwWkR0i8m0RCbnyehH5sVvWCyJyUcb/3+Vuv1NEPp4hv1hEnhWR7SLyMxGpdeUREfmmK98qIqsy9nmP+x87ReRLGfKsx+Ku+6L7/ztE5D0Zct+PZbSI4SH3Wm4TkRU5tlvkXuu97rWPjLS/iBxw67VFRDZmyD8vIkcyrvNbXPlCEYlmyL+esc9vRKRbhjyDYp/XzHWXued7n3tNJMd2d7vbvCQibxppf7c+rRnn+MMZ+2Se/zUZ8m+JyP6MdZeM4rifFpHTIvLzM2zzkYx7akPmtRKRL7nPx+4h9c/1vL3HPdac/1cwqurpB+jOIV8HrDyL8r4F3DKK7VYBP8/4/X+Bf/TomD4AfC3Huj8A03OsuxKYPfScAH8FfN1dvhX4fj7nDPgB8H53+Xrgv7NsMxU4BVRiGiKHgaXuunuBD7nLXwb+wV1eBjzjLl8E7HD3DwFrgSXuuj8C17nLfwH8k7t8B/BNd7kB2OT+9zTgEDDDXfdt4PVnOhbgrcCv3P+uAjYCtWN1LHncG28BngLEvd7P59juceBWd/nrwEdH2h84kO3eAj4PfDKLfCGw4wx1HXY/YZ/XzHUvAFe51+Ip4M1ZtrkA2AqUAYuAV4DgmfYfoT65zv+ozuOQfV4PvD3zvGbZpjZj+SbgaXf5NcDvgaD7eda9Rjmft2zX0evPuJj4ROQ/RGSjq63/MUN+v4jscltUXxGR12BO4pddjX/OKMsXoAZod39fISJ/EJEX3e/zXPkHRORHbstjrwxu2X9QRF4Wkd8CV+f4n6VATFXbsq1X1edU9ViWVTdjXtIATwCvz9Vay8EFwDPu8m/c8oZyC/CUqvZiFERMVV921/0K+JOhZanqHmChiMwEzgeeU9VeVU0CvwXe6e5zHrB+hLJagNPASmAx8LKqtrrbrc22z5BjuQD4raomVbUH81K4cQyPZbTcDHxHDc8BU0RkduYG7rW9HnOtwVz7d4x2//FmMjyv7jmvVdVn1bx5v8PANcrkZuAxVY2p6n5gH3BFHvv7hqo+A3SNsE1nxs8qIB1KSIFyIIJRvmHgBGd+3nzHDwVVIYNNBu/Jss1nVXUlsBy4TkSWi8hUzEvjQlVdDvyzqv4BWAN8SlUvUdVXRvjv14rIFkxr/QbgG658D3Ctql4K3AN8IWOfS4D3AK8C3iMije7N9o+YG/0NmBdfNq4GNo9Qp2zMxbRKcF+YHZgbIRvfdM/j5zKU2FYGbpJ3AjUiMnT/W4FH3eU2ICwDJptbgMaMst4F5sUALADmYXoc14rINBGpxLT00/vswLyIAP50SFk3i0hIRBYBl7nr9gHLxJigQpgHN3OfbMeyFXiziFSKyHTgde4+Y3Uso6X/Wro0u7JMpgGn3Ws9dJsz7a/AL0Vkk4jcNqTM/+Mqhm+ISH2GfJH7Yv+tiLx2FPW3z6thLubcp8l2HdPbZbteI+3/J+71ekJEMu+xclf5PyeuOT6D+9x9vioiZTnqnTcicoeIvAJ8CbgTQFWfxTQQj7mfX6jqbs78vPmOHwoq6t6c6c/3s2zzbhHZDLwIXIi5oTqBPuAREXkX0HsW//079z8bgW9iLgBAHfADEdkBfNX9zzTPqGqHqvYBuzAvtVcD61S1VVXjQLZjAGO+a82x7kxk6y1lC4r4v1X1VcBr3c/7XPknMS+KF4HrgCNA+uWXbg2+CvgFgNuiuxX4qoi8gGllpbe/H6h3XxQfw1yTpHtzfhHTYnoa8/JP7/MXwB0isgnT8o278m9gHsyNwL9izClJVW0HPoo5j7/DmK7SZWU9FlX9JfCkW8ajGJNDcgyPZbSM5lqeaZszrbtaVVcAb3breK0r/w/gHMzL+hjwL678GDDffbH/NfA9GXlMzT6vhtE+k7m2O9P+PwMWuop8LQPWEzDXayXwv4B/zeh13o0xU1+OMdd/Oke980ZV/5+qnuOW+fcAInIuxtIwD6NYrxeRa0d43nxnzE18bsv6k5gxiOXA/weUu63LK4AfYlrYTxf4V2uA9AP9T8BvVPUijI22PGO7WMZyCjNGAdlvzqFE02WJSDCjFXrvCPs147ZC3B5FHWa8aBCqesT97gK+hzk/qOpRVX2X+yL6rCvryNj13cCPVTWRUdazqvpaVb0CY9La68o7VfWDqnoJ8OfADGC/u+6/VHWFql7r1i+9zx5VfaOqXoZRHq+48qSqfsJ96dwMTMnY52eq+mpVvQp4KUOe81hU9T63rDdgXgDpfXw/lly4rc/0dZ6TeS1d5gFHh+zWhjHdhbJsk3N/VU1/twA/ZuD6n1DVlKo6wH9myGOqetJd3uQey9IzHc9ITKLntRlz7tNku46Q+3rl3F9VT6pqut7/ibEs4K5Lb9OEGfe71P19TA0xjPK+YhTHly+PMWCGfCfGDN6tqt2YMbQr3bpkfd7GgvEYg6oFeoAOMeMDbwYQkWqgTlWfBD6OaR2C0dg1Z/E/1zDwsqnDtMzBDFiOxPPAKtckFMaYfrKxGzgXwH1hpFuhWb2EMlgDvN9dvgX4tdtS6cc1k013l8PA2zDmKERkuoikr93dDJhG0ryXAfNeury0B2AZpuX0dff3FHE9yoAPA+vVtVNn7DMfYzp7dIg8gGmBpcuqFJEqd/kNmN7LriH71GOcRB4507G4L5Bp7vJyjHnpl2N4LHNFJD021o/b+kxf56OYa/nnYrgS6NAh447utf0N5lqDufY/dZez7i8iVSJS49alCngjA9c/c4zqnRnyGSISdJcXA0uApqHHkCeT4nl1r1mXiFwpIoJp4Pw0SxlrgFvFeOIuwpzjF860/5DrdZNbD8R4nZa5y9MxJshdmfu4Zb2DgWt8hYh8ZxTnJCsisiTj51sZUDaHMJaMkHsOr8uoZ9bnbUxQj70uMK2aLRmf+135OlyvIIyHym5Ma+xHmJtwNsYLZhuwnQHPrvRFexFj1vgI8JEs/7sKM5azBWPCWc+A58lVwMsYL5V/Ag5oFu8a4OfAKnf5g+4+vwUeJIsXDsYrbCcgOc7FlzAtK8f9/rwrL8d4r+1zj3lxxj5b3O8qjBfcNvc/HmTAW+gWzI31MuZFX5ax/0LMwx0YUpcvu+f8JeDjGfKr3LL2uNeiPmPd79xzvxXX686V3+X+98sYs5pk/PdL7v+sBRZk7POoW9YuXG+2Mx2Le47S2z8HXDLGx7ISY4cf6X4X4P9hXq7byfB8w5go57jLi91rvc+99mVn2t/dfqv72YkZB0qX+9/uttswL8zZrvxP3G23YsZa3j6krusY7sVnn9eB9SsxiuAV4GsZ98JNwL0Z233W3eYlMjz9zrD//824Lr8Blrny17jnbqv7nekd92tXtgP4H6A643lZnaP+v8OYMKOY982bXPm9wE3u8oNuXba4dbnQlQeB1e513gU8MNLzlnEdffPiS59Ay1kiIg8CP1PVteNdF4t3iMj/AQ6p6poRNy4RRGQdxj1940jbTlRK/XkVkS9jpmJsG++6gJnPhrmn3uZL+VZBFYZr9nj1RHqRWSYeIvIbTI/m7aq6dbzrM17Y59U7xHh8/gOwSVXfN9L2Z/UfVkFZLBaLpRixsfgsFovFUpSUtIK68cYbFeNeaj/2U8qfs8Y+A/YzQT5ZKWkF1daWNcKQxTJpsM+AZSJT0grKYrFYLBMXq6AsFovFUpSERt6k9Fi3p4XV65s43N5LY30lt1+7mFXLPEk1Y7FYLJYxYsL1oNbtaeFTT2zlxUPtHO+I8uKhdj71xFbW7WkZ76pZLBaLJQ8mnIL64tN7ONUTJ5ZySDkQSzmc6onzxaf3jHfVLBaLxZIHE87Et6+lm5RmxL5XE2xsX0v3GfayWCwWS7HhWw9KTCKx34jJb79TRO7KWPcxEXnJlWdmxbxbTI77l0TkTWfzv8l0ZAzJ+GTKLRaLxVIS+NmDSgJ/o6qb3bQBm0TkV8BMTNrk5aoaywjlfgEmMdaFwBxgrYgsVdVUPn8aFEgqDNVHwXwSqlssFotl3PGtB6Um4dZmd7kLE659Liaz6v3qJvBSk4wNjNJ6TE3Stf2YtAR5J+la0lBDAEgnRxcxB7mk4WxS1FgsFotlvBgTJwkRWYjJFPk8JsPna0XkeRH5rYhc7m42FzicsVuzK8uLT9+4jHAo0N+DUoVwKMCnb1xWwBFYLBaLZazxXUG5mTd/iEl01YkxK9Zj0gl/CnjczRqZzQg3bOBIRG4TkY0isrG1tXXYDj/d0kws6QySxZIOP93SXPCxWCzFwEjPgMUyUfBVQbmpg38IfFdVf+SKm4EfqeEFTLbZ6a68MWP3ecDRoWWq6sOqulJVV86YMWPYf67ZdjxrXXLJz4aH1r7M8s//gnP+7kmWf/4XPLT2Zc/KtlhGYqRnwGKZKPjpxSfAfwG7VfWBjFU/Aa53t1kKRIA2TOrqW0WkTEQWAUswKaXzIuVk99bLJc+Xh9a+zIO/3kc0kSIUgGgixYO/3meVlMVisXiMnz2oq4H3AdeLyBb38xbgG8BiEdkBPAa83+1N7QQeB3YBTwN35OvBNxY8smE/qkrKUeJJ862qPLJh/3hXzWKxWCYUvrmZq+oGso8rAfxZjn3uA+7zq05e0NWXHDQwlnbG6OpLjkt9LBaLZaIy4UId+U0uQ6GdBmyxWCzeYhWUxWKxWIqSCReLr9SxqUIsFovFYBVUEZFOFdLVlyTpOLR1xfjUE1v58i0Xe6KkrPKzWCylhDXxFRFffHoP7b0JFAgFAyjQ3pvwJFXIuj0t3LNmJy1dfUypCNPS1cc9a3baPFkWi6VosQqqiGhq68FxlFjSoS/hEEs6OI7S1NZTcNmr1zcRDgqVkRAi5jscFFavb/Kg5haLxeI9VkEVEYmUgzNE5rjyQjnc3ktFODhIVhEO0tzeW3DZFovF4gdWQRURuYJdeBEEo7G+kmhi8LznaCLFvPrKwgu3WCwWH7AKapJw+7WLSaSU3ngSVfOdSCm3X7t4vKtmsVgmKOv2tPDeh5/jmi/+mvc+/FzeY97Wi2+SsGpZA/dixqKa23uZZ734LJZJj5+evWnHrHBQBjlm3Quj/g/bg5pEbGs+zc6jHRzt6GPn0Q62NZ8e7ypZLJZxwm/PXi8cs6yCmiQ8tPZlHli7l86+JClH6exL8sDavTYKu8UySfHbs9cLxyyroPIkV/TbXPJi4d9+vTcvucVimdj47dnrhWOWn/mgGkXkNyKyW0R2ishdQ9Z/UkRURKZnyO4WkX0i8pKIvMmvuhVCqQaLTeTwVM8lz5dCB0MtFsvY4rdnrxeOWX72oJLA36jq+Zj07neIyAVglBfwBuBQemN33a3AhcCNwL+LSHBYqZaiw0apsFhKj9uvXUxHNMHeli72HO9kb0sXHdGEZ569q5Y1cO9NF9JQU05HNEFDTTn33nRhXk4YvikoVT2mqpvd5S5gNzDXXf1V4G8Z3PG4GXhMVWOquh/YB1yR7/9WRrLrtFxyS+GsXt9EIpXieEcfL53o4nhHH4lUykapsFiKHAFQUFVQ/4YqztbCNCZjUCKyELgUeF5EbgKOqOrWIZvNBQ5n/G5mQKFllnWbiGwUkY2tra3D/usj1y4mMOQsB8TIJzPBoSdlBHk+7G3poq0rTtJRggEh6ShtXXH2tnQVXLZlOCM9AxbLaFi9vonaijBLZtZw/uw6lsysobYi7FnD0gvLiu8KSkSqgR8CH8eY/T4L3JNt0yyyYYpXVR9W1ZWqunLGjBnDdrjzhqV8/PVLqC0PEQwIteUhPv76Jdx5w9KCjqPUWTKjKi95PsSTDilVEimHWMIhkTK/40mPBrgsgxjpGbBYRoPfThKr1zfRGY2zv62HHUc72d/WQ2c0npcC9HWiroiEMcrpu6r6IxF5FbAI2CoiAPOAzSJyBabH1Jix+zzg6Nn87503LPVNIZWFAsSyvHjLQsXtEHn+7Br2nOjOKi8UR51B4ZhU099WQVksxUpjfSUtXX1URgbUgJdOEtuPtNMdG3gHOAqno0l2HGkfdRm+KSgxGui/gN2q+gCAqm4HGjK2OQCsVNU2EVkDfE9EHgDmAEuAF87mv/2cHZ2rV1DsvYVf7Mrerc4lz4dYDlfAPq9cBC0WD7D50AZz+7WL+eQTWzlyOkrKNc9Xl4X43Fsv8KT8aCL7yFNvDnk2/Gz2Xw28D7heRLa4n7fk2lhVdwKPA7uAp4E7VDWVa/tc+O1RVqpu5r3x7KcylzwfcunmItfZlkmE9TTNjp9OEk6OKNe55NnwrQelqhsY4XhVdeGQ3/cB9xXyv6vXNxFPpjjZnSSecogEA9SUh1i9vsmT1pKQXRkV+0RdPylVpW2ZPGRGTQCojITojSc9ey+UImkniVl1Ff0yL8+JF++FCRcs9uUTnSYrrWkQkEyliCZSJD3IqQQQyTEGFSnyMSg/qYwEs/bErGu/pVg43N7LlIrwIFmp5EPzyzRZCudkwr1Ve+MpHB3Q0ooZnOvxwJQFEAlm7yvlkk8G3nRB9ocll9xiGWtKNR+an6bJxvpKTvbEaGrtZs/xTppauznZE/PsnHgxJ3XCKaika98UQGTA9Jb0IusfICLDzHniyicrxzvjRIbcSZGAkVssxUCp5kPzM6DrVYuncqIzRk88RSKl9MRTnOiMcdXiqR7UHAI5jHm55Nm3nWAERAiI6TmlzXwBMXIvUNVhp1ddeaHkclUvdhf27UfaiQ+xesYd8nIntVj8xIuwO+OBn3OVntx+bJBjRNph4sntxwouG6B76EthBHk2JtwYVENNGYfbo4N6OY4auReIqwDBKMC03vOiB1VdFiSRGjynKCBQU1bcYzk9sRw3Yg65xTIerFrWUPQKaSh+zlXaf7LXNK7d35ohLxaKu2l+FlRFgqSHg9InPChG7gXhoBAMCOFAgLJQgHAgQDAgnoxBLZ1ZS215qF8BBgRqy0MsmVlbcNl+Yr34LJMdv6L5+2maTCSdrNagRBHND5lwCqo7nmJefQWVkaBruw0yr77CMyeJpTNrmVYVIRQUUqqEgsK0qognSuSqxVPp7Ev296Achc6+pGc2YYvF4j1+OjKsWtbALSvm0toVY/fxLlq7YtyyYq43PcESSG434RRUY30loWCAxTOqWdJQw4JpVQQD4mmOk0goyKy6cs6bWcOsunIioaAnLZqndhwf5NyR/n5qx/GCy7ZYLP7gpyPDuj0tPLH5CDNqyjh/Vg0zasp4YvORSTPBeMIpqMwuseM4dPUliCYc/vSyeXT0JgoOSeRni6aprYdgQCgPBSkPBykPBQkGhKa2noLLtlgs/uCnI4Ofyq+6LJTVI7m6rHhcEyacgsr01unsSzKtqoy7rl/CJfOncLInRnN7L4dP9XKyO0Y0nsrb+26yt2gsFstg/Jxj5afyu3B2TdYxqAs9CCANDEt7NJI867ae1KTIWLWsgQ9cvZA/uWweH121mJWL6getT6QcOqIJjnVEOXSql5bOPrpjSVKjmCvlZ4tm0bRKHDWxqlQVx1EcNXKLxVKc+OnI4Kfy29LckZc8X+bUleclz0bx9OU85rEXDvGbl0wyt7qKMCvmT+GyBfVctqCembUDJyjlKN2xJN2xJADl4SBVkRAVkWDW8EV+hgf5zJvP55NPbO1XlsGAMKUszGfefH7BZZcyNgq1pZhZtayBezGN1+b2XuZ5eI/efu1i7lmzk954kopwkGgi5Znyi+bINpBLni+5TIX5mBAnrIJaPKOal453cbSjj45ogt+81NqvsBrrK/qV1SWNU6jKOGF9iRR9iRT0QDgYoDISpDISojwcQER8nZewalkDX7nlYl9u9FJl3Z4W7vr+i3RGkyhwpD3KrmMdPPieSyf1ebEUJ15PrVi1rIFbmk/zyIb99MRTVEWCfPiaRSVx77d2xwgIw+Z1tnXHRl2Gn/mgGoHvALMAB3hYVR8UkS8DbwfiwCvAB1X1tLvP3cCHgBRwp6r+4mz//3Nvu4CPXLeYPce62HiwnU0H23nxcDs9sRSH26Mcbo/yky1HCQicP7uWyxbUs3JBPefPru1Pg25MgcYcGBDjsv7nVy7gC0/tBrxv0UBpTib0k8/9dAcd0WT/bwU6okk+99Md/G7Z9eNXMYvFZd2eFj726GZ63DigR09H2XHkNP/23hUFP8uZY97z3ffNE5uPsHzelKJ/T0QTg4MOgFFWvXn00PzsQSWBv1HVzSJSA2wSkV8BvwLuVtWkiHwRuBv4tIhcANwKXIhJWLhWRJaeTU6oNCLC3PoK5tZXcPMlc0g5ykvHu9h0sJ2NB9vZdayTlKPsPNrJzqOdfOfZg1RFglzSOGAOnFdfgYjgqDEFnje7hjtWncvjmw5zvKOPxqmVfPS6c4r+ZilVmtujeckthofWvjys1e1XlunJzt//eBtdsYHXlKPQFUvx9z/exoa7byiobL/TB/mJFzno/MwHdQw45i53ichuYK6q/jJjs+eAW9zlm4HHVDUG7BeRfcAVwLNe1SkYEC6YU8sFc2p531UL6I0n2Xq4g01uD+vgqV564il+/8pJfv/KScCESEorqxXzpzClMsIVi6dyRcbk2XAwQFt3jKoMU2CxUVseorMvmVVezNgoFfnz0NqXefDX+wgIhALGBP3gr/cBWCXlA0c6jckq87FXHZAXwssnOunsSxJACIqQTCkne+IkU50Fl10Kue3G5O0kIguBS4Hnh6z6C+D77vJcjMJK0+zKhpZ1G3AbwPz588/4v1MrI1RGgvQlHKKJFPGkM8itvDIS4qpzpnHVOdMAaO2K9SurzYfaae9N0NIV46kdx/sny57bUM1K1xx40dw6IqEAiZRDIurQ6ZoCKyJBKiJBKsNBQsHicJT88DWL+l9aabuwo0ZuKS1GegYe2bDfVU7m3gsIJB2HRzbstwrKB3LNVPEgfjSJlJJylERGYQGBeKrwws+bWc2eE91Z5cWC7wpKRKqBHwIfV9XODPlnMWbA76ZFWXYfdhVU9WHgYYCVK1ee8SoFAsYNvDLSvy99CYdYMkVfwqEvkcLJuPAzasq48aJZ3HjRLBxV9rf29I9fbTvSQTzpsK+lm30t3Tz2x8OUhQK8am5d//jVohlVAPTEkvS4XoFlYaOoKsuClIXGL+jrnTcsZX9bN2u2HSeRMh6CNy2fZV9YJchIz0BPPAWqRFMDppSgeJcTzU/8NE365Q2aTtg5VCF5kbAz6U41ycRRb9IHnT+7JquCOt+jeVBe4KuCEpEwRjl9V1V/lCF/P/A24PU60KVpBhozdp8HHPW4Pv29mzTxpENf0njuxRIOCTfzbkCEcxqqOaehmvdc3kg86bDjSEe/wtrb0k0s6bDRHc9aDdRXhlkxv77fJDijpoxYIkUskaK917RoK8uC42IKXLenhU2HOlg4rbLfuWPToQ7W7Wkpelu2JT9CASGWHPwCSymUeZRU068XvZ+myXV7WgZN4WjrjvHJJ7bylVsuLrjub7qggR9vGZ6iwouEndEcjYpc8nz4xc7swQVyyccDP734BPgvYLeqPpAhvxH4NHCdqmZOHloDfE9EHsA4SSwBXvCrfmkioQCRUIDacjO3KeWoUVZJp/9bVYmEAqxYUM+KBWbSb0dvgs2H2vsdLlq6YrT3JnhmTwvPuFElFkytZMWCei5bMIVLGqdQGQnROcQUmHZjD+YzvfosyJxgDMa82RtPlsRgqyU/YjnCeeWS50M6MGo4KIMCo94LBd9Hfpom739qN6d7EwTFjOWoA6d7E9z/1O6C6328M86UilB/oOd0FgIvEnb6OQbb604AHjp21pvwpqc91MU8Uz5a/OxBXQ28D9guIltc2d8BDwFlwK/cHsRzqvoRVd0pIo8DuzCmvzsK8eA7W4IBoaosRJWbPkpV+5VVplmwrjLM65Y18LplDagqze3R/vGrFw+fpjee4uCpXg6e6uXHLx4xDhqza1np9q7Om1WDE1PXFBgbcYJwofg5wdgyeVi9volEarBXWW2FN15lPfEUAZRYMtWfa80r0+T+k71m/NV9O4qAOupJ7qPD7b3UVYRJpLT/nNRVhIv+2RJxk7rqcLkXeDE256cX3wayjys9eYZ97gPu86tOZ4OImMCtGfGwYskUffEB02DKgcaplTROreQdl84lmXLYc3xg/tVu1519+5EOth/p4Jt/OEBVWZBLG03vauWCqcyZUp51gnCFR3ms/JxgXAreQLmwUSryY29LF6e646j7cks6KaJuyvBCKQsG6E2k+u8bVUgoVIa9abA5qiQzlF9AIOhFotFIkH2tPf29s2RKOXK6j3PdMeliZWpFmJO9iaxyL/Ci91fcPsZFSlnIODzUYS7k0HEsgIvm1nHR3Do+8JqFdMeSbD18ur+Hdbg9Sk8sxYZ9bWzY1wbArNpyViyYwsoF9Vw6v566ijAd0QTBgDEFVkVCVEaCZz1u5WfIlPNmVrHnxPCI6+fNLPwBTbfysskLxU9z1USlJ5bEgf63jKpZTDsFFcLUqjC9p1PDXmBTqwp/YQ7NtK0KSYXZdYVn2hYRVJV4+mSICXJajNNNMpEctrZc8vHAKigPGDqOlUgZU2DUVVjVZSGuPnc6V587HYATnX0Z7uyn6YgmON7Zx5Pbj/PkdpMTaunMGi5bYCYMXzinju5Q0jh5uB6BVXmOW/kZL+zQqb685PngZ++slCdBjhexHFEAcsnz4XQ0+5hNRw55PlRFggQwIW3SBPAm03Zrd6xfUQPgLrfmEdJnPGjrzn5ec8nHA6ugfCAcDBAOBqhxFVYy5dDXP46VYmZtOW951Wze8qrZOKq80tLd72yx/UgHiZTy0okuXjrRxfdeOEx5KMDyeXX93oGLplfRJkJZOEiVawoczbiVX2GU/Bxs9XOOiZ+TICcqfg7a98SyK7nuHPJ8aO2OMbQUh/ziwuWiN57KWnY+ERMs2bEKagwIBQNUBwP9UXwdJ8PxIpli6axalsys4dYr5hNLpNie4c7+SmsPfUmHFw6088KBdgCmVkVYMd+YA1csqGd6dVnWwLZjhZ+DrX6+ENPjJpkD546jnkyCnKj4eT189VjLoSy8cMDoy1FGLrll9OStoEQkAFRnTrq15EcgMHg+1lBPwSsWTWPlQhNK6VRPnBcPtbPpoBnDau2OcaonztrdLazdbdzZF06r7O9dXey6s1eOYTSLaZVh2nqGD7ZOq/RmsNUvIqEAPbEkfU5q0NiBH16UlvEl18RWLya8+qlYvXDVLmVGpaBE5HvARzBRxjcBdSLygKp+2c/KTRZyegq641cNNeW8/vyZqCqHT0XdycGn2Hq4g2gixYGTvRw42csPNx8hFBAunFPbr7CWzqzpd7KoiAz+D6/IlZU432zFY82M6jLae+L9YwbiKqkZ1YUPnFuKi1SOXnEueT74qaBy6U8P9GpJMNoe1AWq2iki/xvjJv5pjKKyCson0p6CVAw4XkQTKWoqwpzTUM27Vhh39t3HBqKz7zneSdJRtjZ3sLW5g2/8/gA15SEuzYjO3ji10p0gHKIiHPRkgvDJ3uweXLnkxUKpKlZL/tigw6XJaBVU2A1b9A7ga6qaEJGivbafeGwza7Yd789Ke9PyWXz11hXjXa2CSDte1A5xvKivKmPFgno+kHLo7kvy4uHTbD7YzqZD7TS3R+nqS7J+bxvr9xp39tl15f2ThdPjV35OEC5mWrtjOBneV4ppmXoxcG4pLqyCKk1Gq6BWAweArcB6EVkAFOUY1Cce2zwoLlbKUff35pJXUpkMdbxIh2iaW1/B689vIJ50ON7Zx6YDA9HZO/uSHOvo42fbjvGzbccISNqd3QS7vbhxCnUV4XFxtBgPogln2AtKyS+hmsVi8Y9RKShVfQgToijNQRF5nT9VKow1247nlH/11jGuzBgyEKJpQGHNrC3nvJk1vOuyeUTjSfa1dLPxgDEH7jxq3Nn3HO9iz/Euvvv8IcrDAS6eZ8yBly+cyvmza6gqC40YKzAcgGzvdI8CAPiGFwnVLBaLf4zWSWIm8AVgjqq+2c1+exUmGGxRkcoxephLPlEZqrAcR5kzpYLLF06lL+nQ0ZtgW/NAdIumth76Eg7P7z/F8/tPAa8wrTrCZW509qvOmUZjfWVWR4uPXb+EB9buHVaHj12/ZCwO1WIZkXBASGR5B4QnizvcGKFq0oN49b4drYnvW8A3gc+6v1/GJBosOgVlyc5AbixzyWfXlrNweiVvvHAW0USK4x1RNrrmwE2H2jnZHedkd5xf7jrBL3edAGDx9CouW1DPFYum8urFU6mvjFARDvZHmrYpxgdj4/wVD+c2VLP7eFdWuWU4fQkTs9BRNR8nY1mNIkq5y46jqEJKtd/ByCtHo9EqqOmq+riI3O3+eVJEitIOUhaCbGHBijyz+ZgzNJnj7Npyzp9dyy2XzaM3nmSvG91i08F2thw+TV/Coamth6a2Hn6wqZlwULhwTh0rF9Rz5eKp/K8rF3DbdedQFpr4Y1ejYd2eFj726GZ64ikchaOno+w4cpp/e+8Kq6Ry4GdYq2WzsiuoZbMmhoJK91xyKRTHGaxcks6Zx1kf++Oh/uwN6RijfRmJXvuzO2TEII1mpCnq8yhlx2hf2z0iMg33/hGRK4GOM+0gIo3Ad4BZmMgfD6vqgyIyFdP7WohxvHi3qra7+9wNfAgz3+pOVf1Fvgc0o6qM5o7hXljTq7yb2zIRW8aZPaxplDFnSiWvXjSNaCJFZzTB1sOn+6NbvHyii0RK2XL4NFsOn+aRDfupjASJBAM4qsyfWslfrTqHN144q+BJwqU6UfHvf7yNrtjAQ+oodMVS/P2Pt7Hh7hvGsWbFSyQUyJqzygvv0ie3Zx+bfnL7+I5ND1Us6V5LMuUQHcFZ5z/XNxGNm5if6difacURS6Z/D84gnrn+THx+zS4vD/OsGa2C+mtMQsFzROT3wAzglhH2SQJ/o6qbRaQG2CQivwI+ADyjqveLyGeAzwCfdse1bgUuxCQsXCsiS/PNCXUqmn3uTS55vkyWCNiZY1jTq8uYP7WSN1w4i75EipauPl7Yf4pNrsPFsY4+euMpejGXasfRTj722BYuX1DP2y+Zw9XnTKOhtpyKcP7R2IM5FJRHyWF9I1sj6UxyS+5GhxeNkViOCbm55LlQV4mkHKNARori/vXfvkI0nqTXVSS9cVeZxAcUSN+QnkdagYw0jHPfk7vzqns+VJcZT97ycJDyUJCy9HI4QEXYzNEsz5CVu7KK/uUAf/vD7QXXY7RefJtF5DrgPEyP+yVVHR7bZvA+x4Bj7nKXiOwG5gI3A6vczb4NrMNM/L0ZeExVY8B+EdkHXAE8m88B5epaetXl9DNhWzGT6dY+vbqMZTNreffKRqLxFH/+jed5pcXEDEyTcpTn9p/iuf2nCAgsm1XDZQum8ppzp3H5gqnUVITMROQRyNWItJ7gE49cPYaRehLZUFV6Yyl6XOVwJh5c+7JpYMUHeiLRQUokrVQGeiaxLFMUsnH/U3vyrvtomV4d6Y9AUx4KDCxnKInycNAolHCA8kFKJcinntiWs+yff+waREAYiFOZbluKGGl6vfkGhv7GZwUlIu/KsWqpmwPlR6P5ExFZCFwKPA/MdJUXqnpMRNJv9bnAcxm7NbuyoWXdBtwGMH/+/NH8vafsbemiozdBICAEA0LSUdq64iRSw+3bE5lMk+DhU9FByql/G4GAmHO061gXu4518d/PHaQiHOTixjquWDiNa5ZM54LZNVRGQv1BWy1nZryfAS9wHIfeuENPPEFPzLz0z8QXntxtlEV8sPLoN2/FUyZjQNwEYB6tAgH4ahYPVK9orK+gLGx6FmWuskj3SioipjdSEUkrEbNdRSREpSv/0Lc35iz7mb9eZRaEfqWQVh6QoVDSSmOIwjmTglo4vTiSLY7Ug3r7GdYpMKKCEpFq4IfAx91wSTk3zfEfgwWqDwMPA6xcuXLYer9jV8WTDilVUintz8yZlntBKY5vpSNCD023AfDM31zH802neH7/STYdbOfAyV6iiRTPNZ3iuaZTPPTrvcyoLuOyhcbZ4tolM5gzpcKXmIEThZGegVGWAdA/aK6kkw+6JqyU0pc0vYqeWJJeVwmcib//yXajOAaNiwyYrtI9DxPFP7/n5eH1TWdzmKNi6czqAcURDlLpTqVIK4/KyGBZZSRIZdmAEnnnv/8hZ9nr/9ZMF/XDcaiuyIMxe8EZFZSqfrCQwt3wSD8EvpvR2zohIrPd3tNsoMWVNwONGbvPA44W8v9+4KgzSNmlX8SqhSuodXta+OQTW+mOJUk5Slt3jE8+sZWv3HKxJ0rKL+XnOGnX0sFyVVgwrYoF06p454q59CVSHD4VZcPeNl44cJLNh05zqidOa3eMp3cc5+kdZiD73BnVrFxYz2vOnVZw3caLsqBkHd8oG+fBszsffXGQmSqbN9bZKBCA/3nukA81NlzSOMVVGAEqwqF+hZFWHpnxJSvLglSGQ1RGAv0TzVd9ZV3Osn/5iet8q7f1aC2MUTtfi8hbMQ4M5WmZqt57hu0FM09qt6o+kLFqDfB+4H73+6cZ8u+JyAMYJ4klwAujrd9YEZAAARmspAICJgtJYdz/1G7ae+L98eGSKSWRjHP/U7sLViR+OndURIJZ7fwVGdlK07EEL5gT5oI5tXzg6oX0xpLsONrJH15pY+OBU2xt7iCWdNjX2s2+1m4e++Phguo1nng1KO81a7b61+Z79aKp/T2RtNLo73FEQkZxREJUhNOKI0i1q0AqI0Gu+MIzOcv+yR1X+1bvyUp5SOhLDr8fy0PFo1RHG0ni60Al8DrgEYwH30jK42rgfcB2Edniyv4Oo5geF5EPAYeAPwVQ1Z0i8jiwC+MBeEe+HnxjxVBzoVfmw32tPWS+vxRIKbzS2lNw2avXNxEOSv9E3cpIiN540hPnjopQILuCOoN7cCQUIBKKcM2S6Vx97jRibnSLPx48xR9eOcmmA8ad/UyndueRDhbNqKI8FLTjV6Pk2iXTzQC6O2ieViJpU1a6d1IRCRovzvSy+/vy+3Irke/fftUYHomlUDTHzDP1ZOaZN/PaRtuDeo2qLheRbar6jyLyL4ww/qSqG85Ql9fn2Oc+4L5R1mlcSOZoAeeS51V2Dk2XLURLvhxu7yWWSLK/rQdHTa9vWlXYk7Gz7hxjE7nkQ+nPh1UX5G3L5/DWV80mlnQ4ejrK9f/y25z7vfXfNtBYX8HKhSayxWvOmc6M6rJJEej2bPnP968kIOJ+rAlqMpPKMVk3lzxfxlJB9bnfvSIyBzgFLMrjfyYMfrux+4bj0No9MDPAUWjtTjCvrnDTpNc3elphLZ4x8iz/w+1RDrcf4ccvHiEgcMHs2v7YgZfOr6e2PGwVVgajce23TA5yPZ4e6SdyFZNP8aNVUD8TkSmYBIWbMYrxP/P4nwmD5mgWaJG///ycwBwJBkk6wxV0JOjvy/Dh913Gc00n+eOBdnYd6yTlKDuOdrLjaCfffvYglZEglzROYeXCqVx9zjTjrRUJWQ9BiwVvFIjfjFZB7QFSqvpDN+LDCuAnvtWqiMkVA7HYk7D62vPLpZx9VtpvvHAWN5w/k75kitauGM++ctJEuDjYzsFTvfTGU/zhlZP84ZWTPPTMXhpqykzuq4X1vOac6f5WbhwJkP0lU+TZTyyWYYxWQX1OVX8gItcAbwD+BfgP4NW+1cziKSKSVYt6Yfoaz7xK6QnDC6aFWDCtygS7TaQ40NbD7/e1sfGASdbY3pugpSvGUzuO89SO7HHZJgrl4SC9WRoetufoD34GuZ3sjFZBpe/2twJfV9Wfisjn/amSxQ80h6NFLnmpEgoGqA0GWD5vCq+aW0ffaxx64km2H+nguVfMZOFtRzpGdA45ejraP9em1CK0J1LZjy2X3FIYNp28f4xWQR0RkdXADcAXRaQMazEoKUrB3pyNQlqnImLcqCNBXndeA689dzq9iRTtPXH+eKCdT/5ga8593/rQ71jhJmtcuXAq8+orMqINFLnCGieTq8XiNaNVUO8GbgS+oqqn3QgQn/KvWhaLYVp1hLbueFZ5vqR7V7XlYeZPrTyjgmrvTfDMnhae2WMCnSyYWsllC4zCumT+lP5kjcWosJwcveJccoulWBltNPNeMuY9ZUYqt1j8ZHpVdgU1vSp/BZXJSArlzuvP7U/W2BNPcfBULwdP9fKjF48QDAgXzK5lpauwls2uocINs1PuToAdT/yOR2mxjIayHPm9yvLI7zXh8szOm1JO8+m+rPLJTKkO5Daf7s0qP5JD7hUfXXUuvfEknX0Jtjd39GcXTruzbz/SwfYjHXzzDweoKgtyaWO928OawtwpFb7WbSTsmIilGAjmiO6SS56NCaeg3r2ykQeyhM9/98rGLFtPHnI48VFElqmsRBPZX6u9OeRekR67mlZdxrz6Sq47r4FoPEVrd4wth9r7Fdbh9ig9sRQb9rWxYV8bADNrvcvebJncFGvQ4dHgxdSWCaegntpxnABm4mw6HYaokd95w9Lxrt64Uapmn1SOCuaS+0E4GKCuIkBdRZiGmjIWT6/ixotmE42nOHK6l80HTWbhzYdO0xFNcKJzfLPm5mqMFHu4wmBAsl7XfFrcE42a8hCxnuG5YWvKC391B3Jkq/bqdHvxzplwCqqprYdQUAgGBuycKcehqa3wgKsWSyAgJohqmXl0ZtWVc8HsOt65Yh7RRIq9J7rYfLCdh3+3f/wqWaKTyW9aPosfbxk+tH3T8lnjUJvioCuWGjbxOuDKC6W6LERn3/BIMtVl3qgFL4YVfHMVF5FviEiLiOzIkF0iIs+JyBYR2SgiV2Ssu1tE9onISyLyJr/qZbF4SSQUoK4yzJwpFSyaVsU1S2bw4WsXj2udQq75ZyAt92B5IfjpwX7zJfOIDKljJCjcfMm8gssuZc/7UFBMxHn348V1BKgty+7Mk0ueL7k8bfPxwPVzLtO3MK7pmXwJ+EdVvQS4x/2NGz7pVky+qRuBfxeRszpLi6ZV4qhxqVVVHEdx1MgtFj8JBITqshANNePrkHPO9CrS77B0rykoRl4ouRyw8nDMyskXn96Do8bLqzwcoCwUwFEjL5RIjk6BR50F3/Dzfdbckd0UnUueL9OrIsMUTID8PHB9U1Cquh4T9XyQGKh1l+sYyJh7M/CYqsZUdT+wD7iCs+Azbz6fKZVhJAApVSQAUyrDfObN559NcZYJTC5313zcYIuRz7z5fOqrIpSFA4SDQlk4QH1VxJNnwMnR58glz4emth5SjhJLOvQlHGJJh5Sjnpjnk072+iVyyIuFUn6fdcdTNE6toCoSJBwUqiJBGqdW0JNHCLSxfhI/DnxZRA4DXwHuduVzgcz0qc2ubBgicptrHtzY2to6bP2qZQ185ZaLubSxnlm15VzaWO9ZynTLxOItF83MS14sjOczoG6XTDI+mfJCSCSdYWMW6soLpRicbc6GUn6fNdZXEgoGWDyjmmWzalk8o5pQMMC8+tH3/sa6g/tR4BNuVPR3Y1LC30B2U3DWO0dVHwYeBli5cmXWbVYtayiJCzhRCApky9fohancT8+uP+4f2sE/s7xYGM9noLosRE8siTLgJRuAfqeRgijRyXrVkQDd8eFKtDribfvfa1V61aJ6nt3fnlXuBbdfu5h71uykN56kIhwkmkiRSCm35zFGO9Y9qPczEJHiBwyY8ZqBzIlK8xgw/xUVpTzY6hd+piC5YsGUvOT5cCSHrT2X3AIfvmYRiBAMCJGQ+UbEyAukVKdC9GZRTmeS58O6PS3c9f0Xea7pJM3tUZ5rOsld33+RdW4IrkK4/bpzCQ/RAOGAkXvBqmUN3HvThTTUlNMRTdBQU869N12YV8NprHtQR4HrgHXA9UB6Ru0a4Hsi8gAwB1gCvDDGdRsVpTpL38/GqZ+BaE/nSKiYS54PpXotx5P0XMJHNuynJ56iKhLkw9cs8mSOYYl2oHy9/z/30x10RJODzKkd0SSf++kOfrfs+oLKXr2+ifnTqqjM8CDpjSdZvb7Js953oT153xSUiDwKrAKmi0gz8A/AXwIPikgIk0b+NgBV3SkijwO7gCRwh6oWeQ710iIUFBJZ7HBeuaz6RVNbDxE7r62ouPOGpb5Mei/VaCd+0tweNQuZ50Az5AVwuL2XKRXhQbKKcJDmdn/DiOWDbwpKVd+bY9VlOba/D7jPr/pMdlLZBonOIC8mHFWSydTAmIdAcDK/tSYopZqt2s+en5+9/Mb6Slq6+gb1oKKJVF5ODH5T2v60E4yaHBPkcsnzQXM8LbnkxcKM6ghJx4xDKOY76Rh5oXgRzNLiHblOuxeXI1ewaC+CSPupRCrcQSLVgU+mvBBuv3YxiZTSG0+iar7zdWLwG6ug8sRPJ4m/fG32GyOXPB9KtXVaXRYadm4Fb8KxhHK8+XLJLf4SyTH/LJc8H/75Ha+icshLvTIc4J/f8aqCy/bznfDR687JS54PXjgx+E2Rz6MuPrxMoDeUZ5tOMbuujM5oknjKIRIMUFsR4tmmU9xZcOn+4aeJo60nbsoZEvy3rWf4NciXlJN9GDuX3OIvQR+jVABUloVwSJJylGBAqCz2MBLA8nlTqC0L0h1P4ajpTVZHgiyfN8WT8ot9Sk7xX6EiY3pVhFPd8WHBGwtNoAdm0HJaVRnTqwfMDqpaVIOW2fDTxBFPOoMGydPx5eIeTN4MBgIksyijTIcMy9jh5IjqkPIg2sPq9U3UVYSZXTeQq8srj7VgwJids8kLZfX6JkLBADDgMxYKBjz1tCtmrILKk3T4jrbueH8vZ3p1JK/wHblorK/kwMnuYT2ohdOqPah5aeKoM2gScNokqVq4gqqOBLNm/KyOjG9G3MlKImWuRbrHjDvumJYXgp8ea0saanjpeNewFD9LGmoKLnvn0Y5BEccdhVO9CXYe7Si47FLANhXzxIvwHbm4avFUWrqM4gsIxFMOLV1xrlo81YOalyYBCRBwPfcEBpal8Fu3O0ejIpfc4i+BgBAKuG7l7os+FDDyQmmsryQ6JFGeVx5rn75xGdOqI5QFA4QCUBYMMK06wqdvXFZw2V1Z0mGcST7RsAoqT/z0fHm26RQzqiNEgiaKcyQYYEZ1hGebijv0jp/ecJFQwDSm3da0qlFUXgycZ+s9nUlu8RcToVsIB0wk83AgAIgnkbv9fG5XLWvgy7dczKXz65ldV8Gl8+v5slexD/OUTzSsgsoTPz1fDrf3Mr26bFDvbHp1mSdmiMqwMVsNDfKZlhfC0obqrJ52SxsKN03OqC7rL7Df9JMpL4D02JbIwCdTbhlb/IzcvWpZA5fNr+PAyV52HO3kwMleLptf5/k4jteKY7KHVrNjUGeBX54vfk6c+8h1i/nq2r2DHiBx5YXy5otm8dKJrkEu6yJGXiiqiogQEemPNJBS9SR69tzaMpo7YsNc7efWFq78LPmTjty9en0Tze29zKuv5PZrF3vyrD209mXWbDtOQCAUEhyFNduOs2j6ywVHxVi3p4VPPbGVrr4kScehrSvGp57Y6kkval59BYfbo0MDSTCvviLXLnmxbk8Lq9c3cbi9l0YPz7dX2B5UEXH7tYvpjCbYe6KL3cc62Huii85owhMzxPJ5U6gpC/ZPegyImQDshbvqs02naKgpG5T3paGmzBPTZHc8xdwp5YSCQkqVUFCYO6XcE6eUf37n8qzn5J/fubzgsi1nx6plDTx625X87tPX8+htV3r2snxkw36jnAIBAhJwv428UL749B7auuP0JR2SDvQlHdq6454kWvynmy+irsI0WNPtqLqKEP9080UFl71uTwv3rNlJS1cfUyrCtHT1cc+anZ4EovUK24MqMvoSKeIpx2TRVIe+hDcD9qvXNzGjtpwFPgSGTJsmZ9R47x6f7lUunjFgLuyNJz3JWrtqWQN/+drFwwKfFlML0uINPfHUsPlUAcGThs7LJ7qy5rF6+URXwWWvWtbAg++51Jde5er1TYSD0m+xqYyEPA8WW2gPzSqoIuL+p3YTTTiEA4F+c1Y04XD/U7s9USJ+udk21leyv62brr4B9/ia8hCLphc+BuVFTplcrNvTwhObjzCjpoz5btlPbD7C8nlTCj7foRxzY0o8WW/JUhUx1zfTb8dRIy+UXOEsvQpz6deQgt/BYr0wfdrHpYjYf7LXuFEHBBEhEBACYuSF4qeb7VWLp9LaPdg9vrXbG/d4P51SMluQIuY7HBRWr28quOyKHM4nueQWf/nwNYvcOI4OjjruN57ksSpV/HwngDF9tvcmUMzkYgXaexN5mT59U1Ai8g0RaRGRHUPkHxORl0Rkp4h8KUN+t4jsc9e9ya96TVb8do9vqBnsHt9Q4517vF/jEofbe4cpDK9akH63qi35cecNS7nr+nOpCAdJOuY633X9uZ6kDcnlCeuFh6yf+B0stqmtx523KAhCQEyDO59UOX6a+L4FfA34TlogIq8DbgaWq2pMRBpc+QXArcCFmISFa0Vk6WTLCbV4ehUvHe8i4aQGzUg/b1bhprJVyxq4F3yxZfsdoskvTyM/vSYTyYyoCC6aIbeMPX7lsfrIdYv512f2Dsr8GxBvPGT9xM93glf4mQ9qvYgsHCL+KHC/qsbcbdLuIjcDj7ny/SKyD5MO/lm/6leMDHXXTispL9y1oTTd49OeRuGgDPI0uhcKPhY/x7cCASGoOqjHFBRvoiJYigs/swz7jZ/BYhdNq2Rfaw/iaP+YuqNw7vTRvxfGegxqKfBaEXleRH4rIpe78rnA4Yztml3ZMETkNhHZKCIbW1tbfa7u2PJs0ylm1g52155Z6427NpiX/Xsffo5rvvhr3vvwc565k/ppKvBznMjP8a0Z1ZFh5ryUepPHaiI/A6XKnTcsZdvn38QrX3gL2z7/ppJQTn7jxcTrsfbiCwH1wJXA5cDjIrKY7BOjs1rrVfVh4GGAlStXTiiLvp+mMj97In6bD/30NPKrBVldFiKAG56JgegdXuSxmsjPgGXi4MXE67FWUM3Aj9SEAXhBRBxguitvzNhuHnB0jOs27vhpKlu9volEKsXJ7sGR0r2a8+Cn+bAUI7z7GfXeYikVCn0vjLWJ7yfA9QAishSIAG3AGuBWESkTkUXAEuCFMa7buOOnqWxvSxdtXXGSbrK2pKO0dcXZ21L4ZEI/uWrxVE50xuiJm/GhnniKE52xoo/w7mfUe8vkwi/TfCngp5v5oxgnh/NEpFlEPgR8A1jsup4/BrxfDTuBx4FdwNPAHZPNgw9Ma+OWFXNp7Yqx+3gXrV0xblkx15OeSTzpwBCXTzxK/OcnT24/ZlIvuL8FQF15EeO3C69lclAK4Yj8xE8vvvfmWPVnOba/D7jPr/qUAn5GNggHhWgCnAyPGoBIsLi9yvaf7CUUlEFZblOO48nkZT8pBRdeS/EzFuGIihkb6qiI8PNmXDqzNks4orAn4Yj8JuUoycy5YUCoyBUr+OvCa5kc+O0kVOzYUEdFhJ+RDW6/djGRUJBZdeWcN7OGWXXlRELBojc5pd21HTdhoaPeuWtbLMWO3+GIih2roIoIP29GP+f8+ElNeZigDB6DCoqRWywTnck+lmlNfEWEn5ENoDRNTl2xJPPqh7trd8eS4101i8V3JvtYplVQRcRkvxmz4Wc+KIulFCjFhqVXWAVVZEzmmzEbfvcqLRZL8WLHoCxFTamOnVkslsKxPShL0WN7lRbL5MT2oCwWi8VSlFgFZbFYLJaixCooi8VisRQlVkFZLBaLpSixCspisVgsRYmf6Ta+ISItbmqNoes+KSIqItMzZHeLyD4ReUlE3uRXvSwWi8VSGvjpZv4t4GvAdzKFItIIvAE4lCG7ALgVuBCYA6wVkaWTMSeUZeKwbk8Lq9c3cbi9l0YbFcRiyRvfelCquh44lWXVV4G/xQSnTnMz8JiqxlR1P7APuMKvulksfjPZE81ZLF4wpmNQInITcERVtw5ZNRc4nPG72ZVlK+M2EdkoIhtbW1t9qqnFUhiZub1EzHc4KKxe31Rw2fYZsEwWxkxBiUgl8Fngnmyrs8g0iwxVfVhVV6rqyhkzZnhZRYvFM/zM7WWfActkYSx7UOcAi4CtInIAmAdsFpFZmB5TY8a284CjY1g3i8VTJnuiOYvFC8ZMQanqdlVtUNWFqroQo5RWqOpxYA1wq4iUicgiYAnwwljVzWLxmsmeaM5i8QI/3cwfBZ4FzhORZhH5UK5tVXUn8DiwC3gauMN68FlKGRuF3WIpHN/czFX1vSOsXzjk933AfX7Vx2IZa2wUdoulMGwkCYvFYrEUJVZBWSwWi6UoEdWs3twlgYi0AgfHsQrTgbZx/P+zxdZ7bBmp3m2qeuPZFDzOz8BEvR7FSqnWG87yGShpBTXeiMhGVV053vXIF1vvsaVU6z0SpXpctt5jz9nW3Zr4LBaLxVKUWAVlsVgslqLEKqjCeHi8K3CW2HqPLaVa75Eo1eOy9R57zqrudgzKYrFYLEWJ7UFZLBaLpSixCspisVgsRYlVUGeJiARF5EUR+fl412W0iMgUEXlCRPaIyG4RuWq86zQaROQTIrJTRHaIyKMiUj7edcqFiHxDRFpEZEeGbKqI/EpE9rrf9eNZRy8oxfsf7DPgN17f/1ZBnT13AbvHuxJ58iDwtKouAy6mBOovInOBO4GVqnoREARuHd9anZFvAUMnHH4GeEZVlwDPuL9LnVK8/8E+A37zLTy8/62COgtEZB7wVuCR8a7LaBGRWuBa4L8AVDWuqqfHtVKjJwRUiEgIqKSIc4Wp6nrg1BDxzcC33eVvA+8Yyzp5TSne/2CfgbHA6/vfKqiz41+BvwWcca5HPiwGWoFvuqaZR0SkarwrNRKqegT4CnAIOAZ0qOovx7dWeTNTVY8BuN+lHuL8Xym9+x/sMzBenPX9bxVUnojI24AWVd003nXJkxCwAvgPVb0U6KEETE2uvfpmTDbmOUCViPzZ+NZq8lLC9z/YZ6DksAoqf64GbnLT1j8GXC8i/zO+VRoVzUCzqj7v/n4C87AWOzcA+1W1VVUTwI+A14xznfLlhIjMBnC/W8a5PoVQqvc/2GdgvDjr+98qqDxR1btVdZ6bcPFW4NeqWvStGVU9DhwWkfNc0esxGYyLnUPAlSJSKSKCqXfRD2wPYQ3wfnf5/cBPx7EuBVGq9z/YZ2AcOev737eMupai5GPAd0UkAjQBHxzn+oyIqj4vIk8Am4Ek8CJFHPJFRB4FVgHTRaQZ+AfgfuBxEfkQ5mXzp+NXw0mPfQZ8xOv734Y6slgsFktRYk18FovFYilKrIKyWCwWS1FiFZTFYrFYihKroCwWi8VSlFgFZbFYLJaixCooS05E5F4RuWG862GxjAf2/h9/rJu5JSsiElTV1HjXw2IZD+z9XxzYHtQkREQWuvlwvi0i29z8OJUickBE7hGRDcCfisi3ROQWd5/LReQPIrJVRF4QkRo3J9CXReSPbjm3j/OhWSwjYu//0sEqqMnLecDDqroc6AT+ypX3qeo1qvpYekN31v33gbtU9WJMbLAo8CFMZOXLgcuBvxSRRWN5EBbLWWLv/xLAKqjJy2FV/b27/D/ANe7y97Nsex5wTFX/CKCqnaqaBN4I/LmIbAGeB6YBS3yttcXiDfb+LwFsLL7Jy9DBx/TvnizbSpbt0/KPqeovvKyYxTIG2Pu/BLA9qMnLfBG5yl1+L7DhDNvuAeaIyOUArv09BPwC+KiIhF350lJIAGexYO//ksAqqMnLbuD9IrINmAr8R64NVTUOvAf4NxHZCvwKKMek/N4FbBaRHcBqbK/cUhrY+78EsG7mkxARWQj8XFUvGu+6WCxjjb3/Swfbg7JYLBZLUWJ7UBaLxWIpSmwPymKxWCxFiVVQFovFYilKrIKyWCwWS1FiFZTFYrFYihKroCwWi8VSlPz/nUuX4CLmVjwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sns.FacetGrid(bands_df, col=\"elast_band\")\n",
    "g.map_dataframe(sns.regplot, x=\"price\", y=\"sales\")\n",
    "g.set_titles(col_template=\"Elast. Band {col_name}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "正如我们所看到的，看起来这种分层方案很有用。 对于第一个分层，价格敏感性很高。 随着价格的上涨，销售额下降了很多。 然而，对于第二个分层，随着价格的上涨，销售大致保持不变。 事实上，当我们提高价格时，甚至看起来销量也在上升，但这可能是噪音。\n",
    "\n",
    "现在与使用 ML 预测模型进行的用户分层进行对比："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAADQCAYAAABStPXYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5pElEQVR4nO29fZxcdXX4/z7zuM/JJrubhCRIVhPCo6KB+oAxIn5Ba0FbVPhaS/1SiZaKtMUH6ldKqfwKarFQfr+aFBFtFcT4lFYgJWKMtDwIgQAhAeJGyELCZpNN9nkez++Pe2d3dnfu7szO3J2Z3fN+veY1d87c+5kzM597z+d8PueeI6qKYRiGYVQagXIrYBiGYRi5MANlGIZhVCRmoAzDMIyKxAyUYRiGUZGYgTIMwzAqEjNQhmEYRkViBsowDMOoSOacgRKRlIg8JSLPisgPRaSuiLbuFJGLptjnT0XkkPuZu0RkUzGfOa7t60Tkao/3rhKRP3G3P+x+dlpE1ozb73QRedh9/xkRqXHlN4jIfhHpn0KHa0Rkr4g8LyLnZckvcdt7WkTuF5EWV/46EfmFK98mIsvy+J45dZmsLRH5qvuddovIrSIiOdr9hvu/PCUiL4jI0ax2n8j6zz6VdcydIrIv67g3ufKPur/Df071fcqJ9f/R/i8iERH5tttPd4rIuqz3cvbfcZ/xsax+8JTb/pvc9+5329wlIt8UkWCO3+MpEfmzKb5jnYj8XET2uG3dmPVezrZE5N3j9BoWkQ/maDtn/896v0lEXhGR27Jk4p6PL7jn1pWu3J/+r6pz6gH0Z21/D/irce8HC2jrTuCiKfb5U+C2rNffBz5Rou9yHXB1DnkIeBoIua9PAk4EtgFrcuz3Rvf1wsz3B94KLMn+vXJ8zsnATiAKrAB+CwTddruAFne/rwLXuds/BC51t88B/i2P75lTF6+2gLcD/+3qEgQeBtZN8RmfAe5wtyNA1N1uAH4HHDfVfw6sA/6z3H18iu9p/X90vyuAb7vbbcATOIN2z/47iS6nAR1Zr5vcZwF+BFyc6/fI4zvWAe/O6pe/Bt6Xb1vAAuAIUDfFfiP9P0t2i/t/Zf9/nwC+CwQyv1vWeyXv/3POgxrHr4E3iMg6EfmliHwfeEZEgiLyNRH5jTuCWg8jo4fbROQ5Efk5TqfOGxEJAfVAj/v6D0TkURF5UkS2isgiV36diNwhjlfQkRmluO99SRxvZSvOSZeLc4AdqpoEUNXdqvp8jv3+F/C0qu509zusqil3+xFVPTDFV7oQuFtVY6q6D9gLnIVzUgpQLyICNAGvusecDPzC3f6l28akTKKLV1sK1OAaGiAMvDbFx1wC3OV+XlxVY648yuydaZjr/X+k/6hqF3AUWMPk/deLkf7jttfrboZw+uG0Uvao6qCq/tLdjgM7gClnHbK4CLhPVQen2G+M/iLyFmAR8F/j9vs0cL2qpl2dugrQpWBm64k3Je7J8j7gGVd0FvAlVT0ZuAw4pqpnAmcCnxSRFcCHcE6K04BP4ozU8+GjIvIU8ArOiOY/XPlDwFtV9QzgbuDzWcesBs5z9fpbEQm7neZi4AzgD13dcvEOnNHgVKwCVES2iMgOEfn8lEeMZSmwP+t1J7BUVRM4HfkZnBP7ZOBb7j47gT9ytz8ENIrIwgI/N0POtlT1YRyDdcB9bFHV3V6NiMjrcDzAB7Nky0Xkaff73aSq2ReoG9wL9zdEJDpN3cuK9X/A6T8XikjI/X5vAZZP0X89vyNZF3gAEdmC44n1AZuy3vojt/9sEpHleeiZaW8+8AeMDsryaevi8XrlaHdM/xeRAPCPwOdy7P56nP/zcRG5T0RW5qv/dJiLBqrWPVkeB15mtOM95noB4HgWf+Lu9yjO1NdKYC1wl6qm3AvWg+THD1T1TcBinE6f+eOXAVtEJCM7JeuYn7ueSTdOJ18EvBP4iTuq6gU2e3zeEuBQHnqFgLOBj7nPHxKR9+T5ncAZZY5HRSSMc4KfARyHM91yjfv+1cC7RORJ4F04F61kAZ+ZTc62ROQNONM6y3CM6DkisnaSdi4GNmW8RwBV3a+qpwNvAC7NjO7d77Ea5+K4APjCNHUvF9b/R7kDZ1D1OPBPwP/g9J/J+u8EROT3gEFVfTZbrqrnubpEcbw6cIzzCW7f2gp8Jw89MwOKu4BbVbUjn7ZEZAnOYGLLFM2P7/9/Dtyrqvtz7BsFhlV1DfCvOL+hb8xFAzWkqm9yH59x3WaAgax9BPhM1n4rVDXj6k47u646E7X/gXOiA/wzzvzuacB6nGmpDLGs7RSOMcn384fGteVFJ/ArVe12pwDuBd6cx3HZx2eP2pbhjDjfBKCqv3W/8z24o21VfVVV/9AdNX/JlR0r4DNHmKStDwGPqGq/qvYD9+GsY3nhOcp0L8S7cC6OqOoBdYgB38YZ4VcT1v9H9Umq6l+63/FCYD7wIpP0Xw8m6z/DOIb0Qvf14azp43/F8dryYSPwoqr+U1bbU7X1ERyDnpii7fH6vw34CxH5HfB1nMFKJjijE2dNDeAnwOl56j8t5qKByoctwKfdkRQiskpE6oHtwMXuHP0S4N3TaPtsnGACgHk4o36AS/M4djuOl1MrIo047n4uduOM/KdiC3C6OJFCIRwv5Lk8jsuwGef3iLpTJCuBx3C+08ki0uru915XJ0SkxZ1CAGdUOjICE5E9BXz2ZG29jONZhdz/8F2Zz8/RxolAM04gRUa2TERq3e1mnCmj593XS9xnAT4IPMvsY070f7ff17vb7wWSqvock/TfHG0EgA/jTFFmZA1Z/SQEvB/Y475eknX4BdntevV/EfkKzm911Ti5Z1suY9aVPNqe0P9V9WOqeryqnoAzS/FdVf2i+/ZPGfUG3wW8MFn7xRKaepc5ye3ACcAO90J0COdi9BOcP+cZnD/mV5kDROR64HFVzTXt8FERORtnQNCJE30DThTSD0XkFeARnHlgT1R1h4j8AHgKeAlnkTsX9wH/lqXbh3BGq63Az0XkKVU9T1V7RORm4Dc4I9N7VfXn7jFfBf43UCcincDtqnqdiFyAEwl1raruEpF7cIxaErjCnSZ4VUT+DtguIglX18x3Xgf8g4gozgXnCvfzWsg9Zeipi1dbOPP9mf9JgftV9T/ctsb/T5fgBHpkj8xPAv7RbVeAr6tqZq3me+6FS9z/4VPMPuZE/8cJ8tgiImkco/Rx93M8+292/3ebXwt0Zk27gRMIstldnwziTIV+033vSreNJE50XabdnP1fnFsnvoRj4HY4fwe3qertXm25x52AM7vxq3Ht5dP/J+NGnHPgL4F+YNIw+WKR/PUyqgkR+QnweVV9sdy65IOIfABoV9Vby63LdBHnPpqrVfUDZVZlzmP9f+bxo/+bgZqluK77IlXdXm5d5gIi8lHgb4EnVPXj5dZnrmP9f2bxq/+bgTIMwzAqEguSMAzDMCqSqjZQ559/vuIsgtvDHtX8mDZ2DthjljxyUtUGqru7u9wqGEZZsXPAmM1UtYEyDMMwZi9moAzDMIyKxG7UNYxZyLY9XWzY3sH+nkGWN9exfm0761YXlHzcMMqOeVCGMcvYtqeLz23ayZMv93Dw2BBPvtzD5zbtZNseXysjGEbJMQNlGLOMm+7fw5GBOLFUmlQaYqk0Rwbi3HR/QakODaPs2BSfYcwy9nb1k9KsxG7qpAPf29VfRq0Mo3B886DEKfj2S3Hq1u8Skc9mvfcZcapi7nITgWbk14hT1/55ETnPL90MYzaTzGSHkaxHttwwqgQ/Pagk8NduBuJG4AkReQCn8NiFwOmqGhORNgARORmnLskpOEXCtorIquwicoZhTE1QIKkw3h4Fc+aKN4zKxTcPyi3stsPd7sOpVbIUp1LljZlCWzpa0/5CnLTvMbey516qrxicYZSdlW2NBABxDZKIc6KvbGssp1qGUTAzEiTh1iY5A6d89CrgnSLyqIj8SkTOdHdbCmSXGO50ZYZhFMAXzl9NOBQY8aBUIRwK8IXzV5dXMcMoEN8NlIg04JQIvkpVe3GmFZtxSnB/DrjHLYqWawJiwqS5iFwuIo+LyOOHDh3yUXPDqEymOgd+9lQnsWR6jCyWTPOzpzpnSkXDKAm+Gii3ZPSPgO+p6o9dcSfwY3V4DEgDLa58edbhy4BXx7epqhtVdY2qrmltbR3/tmHMeqY6BzY/fTDncV7yQrl16wucft0WXv8393L6dVu4dauvVb+NOYyfUXwCfAvYrao3Z731U9ya9iKyCogA3cBm4GIRiYrICmAl8Jhf+hnGbCWVzh2t5yUvhFu3vsAtD+5lKJEiFIChRIpbHtxrRsrwBT89qHcAHwfOEZGn3Mf7gTuAdhF5FrgbuNT1pnYB9wDPAfcDV1gEn2FUFrc/tA9VJZVW4knnWVW5/aF95VbNmIX4Fmauqg+Re10J4I89jrkBuMEvnQzDKI6+4eSYheFMIEbfcLIs+hizG0t1ZBhG3nhNEtotwIYfmIEyDMMwKhLLxWcYRsVgZUKMbMxAGYZREWTKhPQNJ0mm03T3xfjcpp187aI3lsRImfGrPmyKzzCMiuCm+/fQM5hAgVAwgAI9g4mSlAnZtqeLazfvoqtvmPm1Ybr6hrl28y6rkVXhmIEyDKMi6OgeIJ1WYsk0w4k0sWSadFrp6B4ouu0N2zsIB4W6SAgR5zkcFDZs7yiB5oZfmIEyDKMiSKTSpMfJ0q68WPb3DFIbDo6R1YaDdPYMFt224R9moAzDqAi8El2UIAEGy5vrGEqMve9/KJFiWXNd8Y0bvmEGyjCMWc/6te0kUspgPImq85xIKevXtpdbtVnNtj1dXLLxEc6+6UEu2fhIwWt+FsVnGMasZ93qNq7HWYvq7BlkmUXxAf5GNmYCU8JBGROYcj3k/RnmQRmGMSd4uvMou149xqvHhtn16jGe7jxabpXKit+RjaUITDEDZRjGrOfWrS9w89YX6R1OkkorvcNJbt764pzOwu53ZGMpAlPMQBmGkTde2Z+95JXCPz/4YkHyuYDfkY2lCEzxsx7UchH5pYjsFpFdIvLZce9fLSIqIi1ZsmtEZK+IPC8i5/mlm2EY06Nak8UmPCLVveSFUmwwQDnwO7KxFIEpfnpQSeCvVfUknPLuV4jIyeAYL+C9wMuZnd33LgZOAc4H/j8RCU5o1TAMo4Ko1iwV69e2c2wowYtdfew52MuLXX0cG0qULLJx3eo2rr/gFNoaazg2lKCtsYbrLziloCAM3wyUqh5Q1R3udh+wG1jqvv0N4POMHXhdCNytqjFV3QfsBc7ySz/DmK3URXKP67zkRnFs2N5BIpXi4LFhnn+tj4PHhkmkUlWRpUIAFFQV1L+p2ul62DOyBiUiJwBnAI+KyAXAK6q6c9xuS4H9Wa87GTVo2W1dLiKPi8jjhw4d8ktlw6hYpjoHPrW2ncC4K01AHPlcJTj+B5lCXggvdvXR3RcnmVaCASGZVrr74rzY1Vd0236yYXsHTbVhVi5q5KQl81i5qJGm2nDJDGspPEvfDZSINAA/Aq7Cmfb7EnBtrl1zyCYYXlXdqKprVHVNa2trKVU1jKpgqnPgynNXcdV7VtJUEyIYEJpqQlz1npVcee6qMmhbGaxsrS9IXgjxZJqUKolUmlgiTSLlvI4nS7TA5RN+B0ls2N5B71Ccfd0DPPtqL/u6B+gdihdkAH29UVdEwjjG6Xuq+mMROQ1YAewUEYBlwA4ROQvHY1qedfgy4FU/9TOM2cqV567yxSBFQwFiOS680VBlBwSftKSRPa/155QXS1rTY9IxqWaeK9tALW+uo6tvmLrIqBkoZZDEM6/00B8b/Q3SCkeHkjz7Sk/ebfhmoMSxQN8CdqvqzQCq+gzQlrXP74A1qtotIpuB74vIzcBxwErgMb/0M4zZjF8ZAry8gkr3FrY8l3tayUteCDGPUMDhUoUI+sT6te1cvWknrxwdIuVOTzZEQ3z5908uSftDidwrT4Me8lz4Oex5B/Bx4BwRecp9vN9rZ1XdBdwDPAfcD1yhqimv/Q3DyI2fUWXVGmY+GM99KfGSF4KXba5wmw34GySR9sjy6yXPhW8elKo+xBTfV1VPGPf6BuAGv3QyjLnAhu0dxJMpDvcniafSRIIBGmtCbNjeUbQXJeQ2RpV+o66fVKvRzgRJLJ5XOyIbjCdL0k+gNL+LJYs1jFnGC6/1OpVpnUExyVSKoUSKZAnqKkU81qAiFb4G5Sd1kWBOT6xUYf1+Tdfu7xlkfm14jKzSamTN3V5lGLOUwXiKtI6OVBVngXqgBNNZkWBuX8lLPhc47+TcxsJLXgh+Ttcub67j8ECMjkP97DnYS8ehfg4PxEoWJFGK+/HMQBnGLCPpzvELIDI6/ZYsQeU/EZkwnSeufK5ysDdOZNyVNBJw5MXiZ0LXt7Uv4LXeGAPxFImUMhBP8VpvjLe1Lyi6bYCAx2Selzz3voZhzCoCIgTE8Zwy03wBceTFoqoTLi/qyovFK1S90kPYn3mlh/i4Wc94moLCqb3w816le585MCYwIhMwce8zB4puG6B//I8yhTwXlf3PG4ZRMG2NUdLuhSfzSKsjLxZxjV/A9cxGtktg/BqiwZwZMBqjlZ2iaSDmcSH2kBeCnwld9x0edAYX7uvM9r7DtgZlGIZP1EeCZJaEMhefoDjyYgkHhWBACAcCREMBwoEAwYCUZA1q1aImmmpCI0YqINBUE2Lloqai2/YTP6P4/CxVn0imc3rDiQqKjzcDZRizjP54imXNtdRFgu76RZBlzbUlCZJYtaiJhfURQkEhpUooKCysj5TEiLytfQG9w8mRrAxphd7hZMnWRKqRdavbuOjNSznUF2P3wT4O9cW46M1LS1OWvQqKe5mBMoxZxvLmOkLBAO2tDaxe3ER7awOhYKAk00Lr17YTCQVZPK+GExc1snheDZFQsCQj+vuePTgmsCPzfN+zB4tuu1rZtqeLTTteobUxykmLG2ltjLJpxysVX8qjVJiBMoxZhp/TQn6O6Du6BwgGhJpQkJpwkJpQkGBA6OgeKLrtasXPKL6GaChnRGZDtHJujzUDZRizjEyhuJaGKEcG4jTXRbjm/NX8XvtC4sl0URF3c31EP9P4GcV3ypLGnGtQp5QggS4wIeBlKnkuKsdUGoZRMt51YitntS/g4LHhEdmBY0Mj28GAE+wQCgQIBZ2gh2BQCAXcRzD32DV7RA9QFwmVLD3OioV17D00gKQVESdEPq3whpbS3DhajfiZcfypzmMFyQvluHk1dB4dzinPFzNQhjEL+fGOV/jyz55lfl2YBXURmusj7nOYBfURmusiI8/NdWGi40bpIhlDNRq1FwwKLx0ZYH5tGFUdCS0v1Yj+i+87ias37aQ/lhzJrj0/GuaL7zup6LarFT8zjg95ZFv3kheK11RhIVOIZqAMYxbS3R9jMJ5iMJ7i1Ryj2PHUR4I0jzFc4RGjtqDeMWzNdRFa6qP0DMZHpp1EhKFEktbGKIf6YiNGLeSGn4eDkvc9UutWt/H1i97Ihu0ddPYMsqyEeeeqmb7h5Ej+w0RKUU2WWaP8ONQfIyCMqZUVEKdv5ouf9aCWA98FFgNpYKOq3iIiXwP+AIgDvwU+oapH3WOuAS4DUsCVqrrFL/0MYzbz/tOWsGR+DfsODXBkIE7PYIKewTiH++P0DMbpGYgznHW/y0A8xUB8iM6eoUladRAYMT4giMBJi5u45zcvO0YtyzMLBQNjphMzx4WCgZzTietWt815g5TN1Zt2TkjOG0umuXrTTh7/v+8tk1b5MZQYW8gRHGM1WICH5qcHlQT+WlV3iEgj8ISIPAA8AFyjqkkRuQm4BviCiJwMXAycglOwcKuIrLKaUIZROD/Z0cm//rqDgXiK2nCQj7xlGVedu3LMPkPxFEdcY5V57hlwDdlAxpA5r7MvkoqT1y87t9+W515jy3OvTdCjqSY0wTNbMM6ILWyIsrAhQk04SCgQGGPAgq4Rm6u5/rr7c+fz85JXEqWoweVnPagDwAF3u09EdgNLVfW/snZ7BLjI3b4QuFtVY8A+EdkLnAU87JeOhjEbuXXrC9zy4F4EJRSAWDLFdx55ifpoiPXrXk86DWlVmmqV1qYo6TSkVEmnlVRaSY+L8nNC1VOuJ+Z4Y5ntI1lGLSNLpEaP7x1O0juc5KUp0ucIMK82M60YHjFqGY9sYX2EloYobU1RFtRHqAkFx0wlZh5eNNWE6B2eODXWVDN3VzmqobbXjPw7InICcAbw6Li3/g/wA3d7KY7BytDpysa3dTlwOcDxxx9falUNo+KZ6hy4/aF9BARCAWedKAgk02m++8hLXH3+6rw+I2OoUmlF1TFgr2txjFhax76fvU86nWYglhprwHIYtcy0Y8r1whQ4OpTg6FCCfVPoFpBsY5blmTVEaamP0NoUpbUhSltjDQsaIkRDAf707a/jtl/+lmBARtZF0gp/dvaKfH/2WceJixrY81p/Tnml4LuBEpEG4EfAVaramyX/Es404PcyohyHTzDwqroR2AiwZs2aSi9aaRglZ6pzYCCeAlWGUqNTKUEprB5UMCAEEcIFpu9LZwyX6oinlkzrGLm6Bi6ZTnNsMEl3f4zDAzF6BhI88NxBnnz5KIm0EhSYVxdBBHoG4mNSIDlragk6mPwm3oAw4o3VR4L0xUZ/gxULalm9pImn9vfQ1hiluS5C2F0zmwtTiictacxpoE4q0X1QpcBXAyUiYRzj9D1V/XGW/FLgA8B7dPSuwU5gedbhy4BX/dTPMGYjoYAQS461WymFaImKCk5W4TUQEAJI/hcWN81eOq3csvUFHn+ph4BAJOgYoqODcS5/5woufccKegbjdPfF6O6L0z0Q58hAzFk3G0xkrZnFOTqYGBnZphUO9zvBIePZd2SIy//tiZHXoYC404phFtZHWFAfpaUhQkuj45G1NkZpa4yyqKmGebVhwsEAgULuOq0wtuzKfXO1l7wc+BnFJ8C3gN2qenOW/HzgC8C7VDV7Ynoz8H0RuRknSGIl8Jhf+hnGbCVXSfbJ5IWQqfAaDsqYCq/XQ1HRd4GA8O3/+d1ItB+MTk1+/7H9fPH9J7NkXi0sGT1GXe8se7oxnYZ4Ks3h/hiH+mN098fo7o9z4327marifTKtHHKPm4pwUEYCPybjcH+MppoQoWCg4ryyQbeMR7ZaqqPyYhkfYp4tzxc/Pah3AB8HnhGRp1zZ3wC3AlHgAfcPe0RVP6Wqu0TkHuA5nKm/KyyCzzAqiw3bO0ikUhzuTxJPpYkEAzTVhkqSSWIgniKAEkumUHUunJNNTYo4Ieu5piFbG6Nkr7bd8PPdnp/7wF+tpbsvRldfzDVsjqd2eCBrDW0gPibIIpFSutxjJuMtX9lKNBQYMWYL6iO0NERY6HpkzppZDa2NURY3RWmIhqcM+CgVmWwd4zNflcqOemXUKiTTlp9RfA+Re13p3kmOuQG4wS+dDGM8k01XGRN5sauPI/1x1L24JdMphtyS4cUSDQYYTKRGLhqqkFCoCxefMnSyiLWVbY2sbPNed8mslw3H046H1Rujq3+Yrj7HQ9vwq8kTt8aSaQ72DnOwd+obpmvCgZEQ/IX1kxeYHE6kCMjUEYxeLKgNc3gwkVNeCkpRJ2vuxlgacx6/pqtmMwOxJGkYucpkSsoPxIrPbrCgPszg0dSEC9iC+uIvmAsbIjnvHVrYMPkUHWTyFgaJhoLMqwvzhraxUW7f/nXHhJLv4Fxcv33ZWXT3jZ1uPNI/ep/ZkYH4mPuChhNpXj06nFf2j7f8/QNZ4fjOutnC+igtje7ambvd2hihNhIiFJBRgyaCeBg1L3k5MANllIRq9ET8THw6W4l5ZAHwkhfC0aHcN58e85AXQkt9hMP98THGT1x5seQyTuCsU7xzZavncem0E9XYP5x0pwuHOdQXG3l0D8T46ZPecWKFZP9oiIYm3CQ92U3AxwYTBAK4Yfky5nkmmZUGqhovltVMtXoi+3sGmT9uOqNUiU9nK36WNx+I5b7S93vIC+FQfyxnaYlC8sKVmkzEY3O9c0PyiYsnTjNOZqBu+qPTHa/MnWocvc9sYvaP/liS/liS/XkYM4AL/t+HRjJ9jMn8UR+mpT7KwoYoCxrCRIJBNywfx4CJY8gCJSrkNOsMVLVeLKuZavVE/CxlMFvx00D52bZXep1C7g2rND565vIxrzNrZcmUkkil6XM9M8cjG3bvNxvN/vFwx2HPtl86PJhX9o+mWteAjc/+UVeadayCDZSIBICG7JtuK4lqvVhWM9XqifhZysCoLJK54p0nkVcjo2tlzuv5dRGWLxg72FI3HD+ZVlZ/+X7Ptv7i3W8YiV48Mi6lVTIr+8exoQTH8sj+MV3yMlAi8n3gUzhZxp8A5onIzar6NZ/0mjbVerGsZpY317Gvu5++4dHQ48aaECtaKidlihcCoM6Ji0pF5SEzSkfKI8rQSz5bEXHLoUyRIeSz564cMWSptJJMpUmlHc/s6GCCQ/1OSH4mhVV2OP6RwQQ9A04gSLHk60GdrKq9IvIxnDDxL+AYqoozUDMxbWNrXGN5W/sCHnGnCxRIppw6RP/7rMrOlbhhe8dIiG5KR8N1zdueffg5fTgbCQcDnmmuljY7z+lsA5bOGLDR12ff9Mui9cjXQIXdtEUfBG5T1YSIVOR/u35tO5/9wZP0Dg2gZOZJSzdtY2tcE7n3mQM5C5Pd+8wBrjx3VfkUm4IXXuvl6FCCdHrUsA4nUySnSjmQJzaQqRzMQJWeQECIjET1FZi0Md/PyHO/DcDvgHpgu4i8DqjINainO4/SOzT2nozeoSRPdx4tSfsbtncQT6Y4eGyY51/r4+CxYeLJFBu2T36z3mxm3+FBggGhJhykNhykJuxE9uybYpG13Awl0qTSoxcpBVLpwgqqebFtTxef27STJ1/u4eCxIZ58uYfPbdrJtj2Vk+fMMCqdvAyUqt6qqktV9f3q8BLwbp91mxa3P7TPSeHhvlac1B23P1SaZbwXXuvl8ECcZEoJipBMKYcH4rz4WkXa6xkjlXZS1AwnUsSSqZEyCpVMKQqqeXHT/Xs43B8nlkqTTEMsleZwf5yb7t9TdNuGMVfIy0CJyCIR+ZaI3Oe+Phm41FfNpknfcDJnmeG+HMXKpkMmpUvATcmfyWYcn2OLrdm0NkRIufV1FOc5pY58rvJiVx9pRvOOqULalRszT9jjBlMvuVEZ5DvFdyewBSfLOMALwFU+6FM0fs81R0IB0mllOJFiKOF4DOm0EgmV6M60KqSxJjwh+k1c+VzFa7wyh8cxZWV8eqKp5IXgZeLM9BVPvlfVFlW9B9w0XKpJnJDzOUdrQ3TEU4BRj6G1YfLEjrOZQ/0xAjJ6QgpOkEQ+ZQtmLbYqXzB+XuhXL85tiLzkhWB/tX/ka6AGRGQh7m8uIm8Fjk12gIgsF5FfishuEdklIp915QtE5AERedF9bs465hoR2Ssiz4vIedP5Qn6PZvqGEzlTpvQNT8wKPFeIJ9MoY422uvJScOvWFzj9ui28/m/u5fTrtnDr1hdK0q5RWXjNQpRiduLeZw4WJC8E86ByE/KYPvWS59w3z/3+Cqeg4OtF5L+BVuCiKY5JAn+tqjtEpBF4QkQeAP4U+IWq3igiXwS+CHzBXde6GDgFZypxq4isKrQmlN+jmdd6c3sFXvJCuXXrC9z+0D4G4inqI0H+7OwVFR2qDRBPpnKu+8WTxTvZt259gZu3vjjyunc4OfK6kn8XG1UXjtd1qxTLRDGPuVUveSHYf52btEfhJy95LvIyUK6ReRdwIs7A4HlVndRlUNUDwAF3u09EdgNLgQuBde5u3wG24dz4eyFwt6rGgH0ishc4C3g4728zAyQ9flwveSHcuvUFbnlwLwGBUMC5wfiWB/cCpbkY+3Vfzvjy4lPJC+EbWcZpvLySDZRROMMe4f1ecqOy8QrkLSTAd1IDJSJ/6PHWKhFBVX+cz4eIyAnAGcCjwCLXeKGqB0Qkc4VcCjySdVinKxvf1uXA5QDHHz/zmQpKUSXSi9sf2oeqkn0+BtwQ+WIvxtv2dLk3MCdR4JWeIZ47cIxbPnpG0UaqWhOIVivlPgf8wv5rYzxTeVB/MMl7CkxpoESkAfgRcJWbLslzV4/PGCtQ3QhsBFizZs2M993JKnMWS99wckLbpQqR//LPnuXYUBJhVNdjQ0m+/LNn+fXqc4pu35g5yn0OGMZMMamBUtVPFNO4mx7pR8D3sryt10Rkies9LQEyt9Z3Atn545cB3sVQykS1eguZombj28qn2JlROH4OZAyjFPjdRwO4Yd855PmSd7kNEfl9nACGmoxMVa+fZH8BvgXsVtWbs97ajHOT743u88+y5N8XkZtxgiRWAo/lq58xOTZ9MrPY721UOgvqwhwenBhKUKpaTqU4B/Itt/FNoA4nvdHtOBF8UxmPdwAfB54Rkadc2d/gGKZ7ROQy4GXgwwCquktE7gGew4kAvKLQCD7DMKoX8zpnltpoiKZUit6sisVN0QB10dLUsZ0xAwW8XVVPF5GnVfXvROQfmWL9SVUfwrtvvcfjmBuAG/LUyTCMGcZPI2Je58yyvLmOrlCA12WVJhqMJ2lrrJnkqPypCwcZTEz0Meq86njkIN/pwGH3eVBEjsPxcFbk/SmGYcwKljXXAowE28g4uVE9rF/bTiKlDMaTqDrPiZSyfm17Sdo/75Tc0cFe8lzka6D+Q0Tm4xQo3AHsA+7K+1MMw5gV/P2FpzKv1hlxZzybebUh/v7CU8unVJlpjOa+jHrJK4V1q9u4/oJTaGus4dhQgrbGGq6/4JSS1Sw72Btnfm1o5EbrgMD82hAHe/OvtJvvFN8eIKWqP3IzPrwZ+Glh6hqGMRM0RoP0xSZOrTRGiy8qt251G7d89Aw2bO+gs2eQZVVSiLGlIUJ3/8QLY0sJMu5/8p2vH5PtJFte6axb3ebbf7ff7R/ZtxapKp09+deJy9dAfVlVfygiZwPvBf4R+Bfg9wrQ1zCMGeDUpfPZc/AYx4ac0jMBcbyc1YvnlaR9vy5q46syZ8uLpaU+wpGB+ISqzy31xRuohzuOsLgpSt9wkngqTSQYoLEmxMMdR7iyyLYjQclZyicSrPzQkeXNdXT1DVOXtcY1lEixrLku7zby9UEzw7HfB76pqj8D5m6xH8OoYNavbaexJsKKlnpOPa6JFS31NNZESra24Bd+Zmnpj6dY3lxLfSRIOCjUR4Isb65loATFKff3DNLSEKW9tYHVi5tob22gpSFakKfgRVNt7lI2TbWVX8qmFGtc+XpQr4jIBuBc4CYRiVLY/VYVwR/9y/8QCQqRUJBIKEB05DH6uiYcHJWHg9SER+U1ocmnSFSVSTJlGMaMsG51G9dD1U3DBQKSsxJzoAQuVGY03946Wl6jVBFrpfAUvFjZ1kg42E/v0Kh31lQb4oSFxZcJAf/yc0Jp+mG+BuojwPnA11X1qJsB4nOFq1xennipx9f2V1xzL5FQgEjQMWqRzCPrdTTkGsBw1rb73mT8YvdrWQY06B4/sb2gVQg18HdtwS9qw0J/bKKBqgsX36fXr23n2s27GIwnqQ0HGUqkShaxNhNtL54XKnnb2/Z0ce3mXYSDwvzaMF19w1y7eRfXQ0mNVDFt5ZvNfJCs+56yM5VXEx/7veOJJ9PEU2kSqTSxhPuckSWd7UTKeR1Ppkf2jyfTeWXhzRxT6lp9l33n8bz2C4qMGMZoKEDYNY6T8cUfPT3iMWaMXU1onCfpGsVRD3PU6zQqDz9Hxn5x2tJmdh84Ru/w6NpZU02Ik5YUv3bmp1dZrW1v2N5BOCgjnl9dJMRgPMmG7R0V01dKc8twlXDDh07LKVdV1K2Sq6qjhfcUFB2ZA0+m0px63X95tn/7n7xlxNgNJxxDFUukiLmGL5ZMEU8q8VSKWCJ7vxTxVJpHOo54th0JBfIqAJhSZcgtR58vd/9mf977Fspb/59fjHiRkVC2JznRo3SMabZHOLnxe7rzqONVZjzKcJAa1ygHRMjMts61adeZGBn7QcZbWNgQLbm3AP56ldXY9v6eQeaPW8uqDQdLsnZWKmadgZpOgkLJuphNek/8FHdAn3vy4smVm4ITvvhzz/de+Mr7nFIcKSWWdIxePJlmOJEaMXTDiRSxRIph1zAOJ9Ij+37l57s9237vSYsYTqYcg5rlNcbctjPGNJ5MkyykmAtwsHd46p2myQW3/XdOeUAgHAyMMYyRYIDwOEM5Gf9w3+4RLzESFNegZnmV7nYkM1UbzGw7r8tJNYyMc1Gta2fVip9rZ6Vi1hko8cjFMhsG0SJCJORM4TUWeOxkBupfL12TdzvptI4aR9eovfOrv/Tc/x/+8LQRL3IoMdawZjzKWJZxzEyxDrvTr/u6Bwr6nuCEKmfami4bftUx7WPLTTWMjL2oxrWzasXPtbNSMesMlFcF5xJUdjZwIqpqI0FqI0HmMXWo6yVnFVdQbzKv8jdfOpeheNLxHpPOtGksmWIo4U6tut5fLJli2PUoRz1D5d8fecmz7dWLG0e8xkRq7FpkosI7UzWMjI3yUw0eq28GSkTuAD4AdKnqqa7sTcA3cUp2JIE/V9XH3PeuAS7DuefqSlXd4pduxuygtTEKRKd9/GQG6t4r30lalXTWOmTmdSqVJpbKTKk6Bm8o6awlZozlFd97ctp6FUs1jIyNyqDSPVY/Pag7gduA72bJvgr8nareJyLvd1+vc9MnXYxTb+o4YKuIrJpOuY1wgDEl07PlhpEvgYAQKCJH9xWUz0BVw8jYMPLBNwOlqttF5ITxYqDJ3Z7HaMXcC4G7VTUG7BORvcBZwMOFfu5nzlmZMy/WZ85ZWWhThlG1VPrI2DDyYabXoK4CtojI13EC697uypcCj2Tt1+nKJiAilwOXAxx//MT1jSvPXQXA7Q/tYyCeoj4S5M/OXjEiN4xqZ6pzwDBmCzNtoD4N/KWbFf0jOCXhzyV3bHfOlWhV3QhsBFizZk3Ofa48d5UZpFlC0CP9TSkyZoQEkjl6UKjCIz7zOQcMYzYw0yszlzKakeKHONN44HhMy7P2W8bo9J8xh1nZmjvyzEteCLmM02RywzBmlpk2UK8C73K3zwEyi0WbgYtFJCoiK4CVwGMzrJtRgfy2O/e9O15ywzBmD36Gmd8FrANaRKQT+Fvgk8AtIhLCKSN/OYCq7hKRe4DncMLPr5hOBJ8x+/C656jS70UyDKN4/Iziu8Tjrbd47H8DcINf+hiGYRjVhd0dZMxZvOIsrGKJYVQGZqAKxOvaVYprmtefMZf/pIV1udMpeckLIRzM/ct6yQ3DmFnsTCyQhQ25K917yQvBK6HtbEh0O13iqdwJX73khZBK527DS24YxsxiBqpAWuojE360gCsvFq9KFgVWuJhVDObKWzWJvBCCgdzd30tuGMbMYmdigfTHUyxfUEt9JEg4KNRHgixfUMtA3IIO/SDXTbqTyQuhIZK7bpOX3DCMmWXWldvwm0wpg/bWhhHZYDxJW2NN0W17lLIqyfqWMZF+j0GFl9wwjJnFPKgCWb+2nURKGYwnUXWeS1XKoMajYq+XvFLwSjtUinREfuJV0LCYQoeGYZQOM1AFsm51G9dfcAptjTUcG0rQ1ljD9RecUpLM0XXRIEFxgiIE5zkojrzotj2MnJe8EFa1NRQkrxQywScio49suWEY5cWm+KaBX6UMVrY18rvD/fQOJYmn0kSCAZpqQ5ywsPgL/RuXNfHwvp6c8mJZvbiB3Qf7csqLJRQQkjnWm0Il8M6WNkXpPBZDdaLcMIzyYx5UBbF+bTvJlJJKK6rOc7JE04c7O3sLkhfCll1dBckLYV5taMIanADza4sfW33kzONztv2RM62EhWFUAmagKozhRIp4Kk0y7dzrM5wozYL9oNvO+OmswRK079VGKdpetaiJRU3RMVGTi5qirFxUvOf3cMcRmuvCI5kjAgLNdWEe7jhSdNuGYRSPGagK4sb7djOUSBMOBKgJBQgHAgwl0tx43+6i267Wm4DXr20nEgqyeF4NJy5qZPG8GiKhYEm8yhe7+ugbThIOBqgJBwgHA/QNJ3mxa+J0ZaFYVhDDKB47XyqIfYcHCQgEAoKIEAgIAXHkxbJ0nhMGrzr6yJZXKn4GpcSTaRAIiCAIATc6JV6CKD6vu7Tm8D3XhlEwfpbbuAP4ANClqqdmyT8D/AVOWY2fq+rnXfk1wGVACrhSVbf4pdtc5CsfPI0r79pBfzxFWp3prIZIkK988LSi2w5I7mwXpYoy9ysoJRwUhhKQTisio0Y7EixecTNQhlE8fkbx3QncBnw3IxCRdwMXAqerakxE2lz5ycDFwCnAccBWEVk112pCtbfU8/zBPhLpFKruWpHCiSWIhlu3uo1bL3kzG7Z30NkzyLLmOtavbS/Jhb8hGqJ3OJlTXgq27eliw/YO9vcMsryEeq9a1MS+7n76hkejJhtrwqxoqezweMOYK/hZD2q7iJwwTvxp4EZVjbn7ZMK8LgTuduX7RGQvTjn4h/3SrxJ536mLef61vpGRfMZIve/UxSVp3y9P5JTj5vH8wV6ODiVGvLP5tWFOXFx8IMO2PV1cu3kX4aAwvzZMV98w127exfVQ9HdZv7adazfvYvG8ELXhIEOJVMluurasIIZRPDO9BrUKeKeIPCoivxKRM135UmB/1n6drmwCInK5iDwuIo8fOnTIZ3Vnloc7juSMWCtVVNm2PV1csvERzr7pQS7Z+Ajb9hQfBg7Ohb6hJsyKlnpOPa6JFS31NNSES3Kh37C9g3BQqIuEEHGew0Fhw/aOotv2c31rgUc5EC95Iczmc8AwspnpG3VDQDPwVuBM4B4RaSf3wDLndL2qbgQ2AqxZs2ZWTenv7xlkYX2UlobRwAVVpbOn+CAJPz2RdavbuB58mT7c3zPI/NqxF/XacLAkvwn451W2NkY5PJjIKS+W2XwOGEY2M22gOoEfq6oCj4lIGmhx5cuz9lsGvDrDupWdTCLausjo3zKUSLGsua7otjds7yCRSnG4f2yWig3bO0pygfbrQr+8uc637Bp+0h9P8boFtXT3x0f0bmmIWNZ7wyiAmZ7i+ylwDoCIrAIiQDewGbhYRKIisgJYCTw2w7qVHT8T0b7Y1Ud3X5xkWgm66YO6++IluefHT97WvoADx2IMxJ31oYF4igPHYrytfUG5VZuU5c11hIIB2lsbWL24ifbWBkLBQEkGG4YxV/DNQInIXThBDieKSKeIXAbcAbSLyLPA3cCl6rALuAd4DrgfuGKuRfCB44Vc9OalHOqLsftgH4f6Ylz05qUVf8+Pn9zzm5cLklcKfg42DGOu4GcU3yUeb/2xx/43ADf4pU81sG1PF5t2vEJrY5Tj3aiyTTte4fRl84s2Un7e8+Mnr/TGgLEZL1RH5ZWKn+tyhjFXsGzmFUR2xBpAXSTEYDxZknWiar3nJzvkPpe8kvFrXc4w5gqW6qiC2N8zSO24+kyliljzM6edn0RDubuol9wwjNmDneUVxPLmOobGZQAvVRSfn/f8+ElrQ6QguWEYsweb4qsgMpkNBuPJkmc2gCqdchKhtSHM4YHRLBUL68NIpadhNwyjaMxAVRC2sD6RzL1hi+eNepGD8SRtjZWdhd0wjOIxA1VhVKWX4yN+e5WGYVQutgZlVDTVunZmGEbxmAdlVDzmVRrG3MQ8KMMwDKMiMQNlGIZhVCRmoAzDMIyKxAyUYRiGUZGYgTIMwzAqEj/LbdwhIl1uaY3x710tIioiLVmya0Rkr4g8LyLn+aWXYRiGUR34GWZ+J3Ab8N1soYgsB94LvJwlOxm4GDgFOA7YKiKr5mJNKGP2sG1PFxu2d7C/Z5DllhXEMArGNw9KVbcDR3K89Q3g80B2wYQLgbtVNaaq+4C9wFl+6WYYfrNtTxfXbt5FV98w82vDdPUNc+3mXWzb01Vu1QyjapjRNSgRuQB4RVV3jntrKbA/63WnK8vVxuUi8riIPH7o0CGfNDWM4siu7SXiPIeDwobtHUW3beeAMVeYMQMlInXAl4Brc72dQ5azJJ2qblTVNaq6prW1tZQqGkbJ8LO2l50DxlxhJj2o1wMrgJ0i8jtgGbBDRBbjeEzLs/ZdBrw6g7oZRknxs7aXYcwVZsxAqeozqtqmqieo6gk4RunNqnoQ2AxcLCJREVkBrAQemyndDKPUrF/bTiKlDMaTqDrPloXdMArDzzDzu4CHgRNFpFNELvPaV1V3AfcAzwH3A1dYBJ9RzVgWdsMoHt/CzFX1kineP2Hc6xuAG/zSxzBmGsvCbhjFYZkkDMMwjIrEDJRhGIZRkYhqzmjuqkBEDgEvlVGFFqC7jJ8/XUzvmWUqvbtV9fzpNFzmc2C2/h+VSrXqDdM8B6raQJUbEXlcVdeUW49CMb1nlmrVeyqq9XuZ3jPPdHW3KT7DMAyjIjEDZRiGYVQkZqCKY2O5FZgmpvfMUq16T0W1fi/Te+aZlu62BmUYhmFUJOZBGYZhGBWJGSjDMAyjIjEDNU1EJCgiT4rIf5Zbl3wRkfkisklE9ojIbhF5W7l1ygcR+UsR2SUiz4rIXSJSU26dvBCRO0SkS0SezZItEJEHRORF97m5nDqWgmrs/2DngN+Uuv+bgZo+nwV2l1uJArkFuF9VVwNvpAr0F5GlwJXAGlU9FQgCF5dXq0m5Exh/w+EXgV+o6krgF+7raqca+z/YOeA3d1LC/m8GahqIyDLg94Hby61LvohIE7AW+BaAqsZV9WhZlcqfEFArIiGgjgquFaaq24Ej48QXAt9xt78DfHAmdSo11dj/wc6BmaDU/d8M1PT4J+DzQLrMehRCO3AI+LY7NXO7iNSXW6mpUNVXgK8DLwMHgGOq+l/l1apgFqnqAQD3udpTnP8T1df/wc6BcjHt/m8GqkBE5ANAl6o+UW5dCiQEvBn4F1U9AxigCqaa3PnqC3GqMR8H1IvIH5dXq7lLFfd/sHOg6jADVTjvAC5wy9bfDZwjIv9eXpXyohPoVNVH3debcE7WSudcYJ+qHlLVBPBj4O1l1qlQXhORJQDuc1eZ9SmGau3/YOdAuZh2/zcDVSCqeo2qLnMLLl4MPKiqFT+aUdWDwH4ROdEVvQengnGl8zLwVhGpExHB0bviF7bHsRm41N2+FPhZGXUpimrt/2DnQBmZdv/3raKuUZF8BvieiESADuATZdZnSlT1URHZBOwAksCTVHDKFxG5C1gHtIhIJ/C3wI3APSJyGc7F5sPl03DOY+eAj5S6/1uqI8MwDKMisSk+wzAMoyIxA2UYhmFUJGagDMMwjIrEDJRhGIZRkZiBMgzDMCoSM1CGJyJyvYicW249DKMcWP8vPxZmbuRERIKqmiq3HoZRDqz/VwbmQc1BROQEtx7Od0Tkabc+Tp2I/E5ErhWRh4APi8idInKRe8yZIvI/IrJTRB4TkUa3JtDXROQ3bjvry/zVDGNKrP9XD2ag5i4nAhtV9XSgF/hzVz6sqmer6t2ZHd277n8AfFZV34iTG2wIuAwns/KZwJnAJ0VkxUx+CcOYJtb/qwAzUHOX/ar63+72vwNnu9s/yLHvicABVf0NgKr2qmoS+F/An4jIU8CjwEJgpa9aG0ZpsP5fBVguvrnL+MXHzOuBHPtKjv0z8s+o6pZSKmYYM4D1/yrAPKi5y/Ei8jZ3+xLgoUn23QMcJyJnArjz7yFgC/BpEQm78lXVUADOMLD+XxWYgZq77AYuFZGngQXAv3jtqKpx4KPAP4vITuABoAan5PdzwA4ReRbYgHnlRnVg/b8KsDDzOYiInAD8p6qeWm5dDGOmsf5fPZgHZRiGYVQk5kEZhmEYFYl5UIZhGEZFYgbKMAzDqEjMQBmGYRgViRkowzAMoyIxA2UYhmFUJP8/fUvgHUgOSzUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g = sns.FacetGrid(bands_df, col=\"pred_band\")\n",
    "g.map_dataframe(sns.regplot, x=\"price\", y=\"sales\")\n",
    "g.set_titles(col_template=\"Pred. Band {col_name}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我真的很喜欢这个图，因为它传达了一个非常重要的观点。如您所见，预测模型分层在 y 轴上分割个体单元。如第一幅图所示，在第一个分层所对应的日子里，我们没有售出很多冰淇淋，但在第二个分层对应那些日子里，我们确实卖得更多。我觉得这很神奇，因为预测模型正在做它应该做的事情：它预测销售。它可以区分冰淇淋销售量低和高的日子。\n",
    " \n",
    "唯一的问题是预测在这里并不是特别有用。最终，我们想知道什么时候可以提高价格，什么时候不能。但是，一旦我们查看预测模型分层中线的斜率，我们就会发现它们并没有太大变化。换句话说，由预测模型定义的两个分层对价格上涨的反应大致相同。这并没有让我们深入了解哪些日子可以提高价格，因为看起来价格根本不会影响销售。\n",
    " \n",
    "## 关键思想\n",
    " \n",
    "我们最终正式确定了条件平均干预效果的概念以及它如何对个性化有用。即，如果我们能了解每个单位对一种干预的反应，即如果我们能了解干预效果的异质性，我们就可以根据单位的个体特点给予最好的干预。\n",
    " \n",
    "我们还将这个目标与预测模型的目标进行了对比。也就是说，我们正在重新考虑估计任务，从预测原始格式的 \\\\(Y\\\\) 到预测 \\\\(Y\\\\) 如何随 \\\\(T\\\\)、\\\\(\\frac{\\delta y}{\\delta t}\\\\)。\n",
    " \n",
    "可悲的是，如何为此构建模型并不明显。由于我们不能直接观察弹性，因此很难建立一个预测它的模型。但是线性回归拯救了我们。通过使用适合预测 \\\\(Y\\\\) 的回归模型，我们找到了一种方法来预测 \\\\(\\frac{\\delta y}{\\delta t}\\\\)。我们还必须包括干预和特征的交互项。这使得我们对每个客户的弹性预测都不同。换句话说，我们现在正在估计 \\\\(E[T'(t) | X]\\\\)。然后使用这些弹性预测将我们的单位分组为或多或少对干预敏感，最终帮助我们确定每组的干预水平。\n",
    " \n",
    "![img](./data/img/causal-model/economists.png)\n",
    " \n",
    "所有这一切引发的一个自然问题是，我们是否可以用通用机器学习模型替换线性回归并使用它来预测弹性。答案是肯定的，但有一些需要注意的地方。本章使用了一个非常简单的 CATE 模型，因为我认为使用线性回归更容易理解它们背后的概念。不过不用担心。我们将在接下来的章节中看到一些更复杂的模型。但在此之前，我首先需要介绍一个非常重要的话题，即我们如何比较两个 CATE 模型并决定哪个更好。"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
